\u;>  ■xv  - ■.  ■> 


■  V  "  V  ■  '.  »  '.  < 


’*  ■▼  ^  »,■%■  >w-j  i  -v  ’  v  ■■.  ^ v\  » .,,  w ."  w."  y*w^u  ■  «  -  v  - 


m flLt 


*fr 

CD 

CO 


> 

^  *  K. 

v. 

\ 


i 


00 

< 

I 

O 

< 


AFGL-TR-87-0229 


Lg  WAVE  EXCITATION  AND  PROPAGATION  IN  THE  PRESENCE  OF  ONE-,  TWO-, 
and  THREE-DIMENSIONAL  HETEROGENEITIES 


R.  D.  HERRMANN 


i 


SAINT  LOUIS  UNIVERSITY 
221  NORTH  GRAND  BOULEVARD 
ST.  LOUIS,  MO  63103 


15  July  1987 

Final  Report 

1  April  1985  -  31  -  March  1987 


DT!C 

ELECTED 
NOV  2  5  1987 i 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


AIR  FORCE  GEOPHYSICS  LABORATORY 
AIR  FORCE  SYSTEMS  COMMAND 
UNITED  STATES  AIR  FORCE 

HANSCOM  AIR  FORCE  BASE,  MASSACHUSETTS  01731 


**  U  1 


o 


'»  'v.-T*v 


062 


E59S53 


SC C'J*'T>  ClaSSi F 'CATiON  Of  TNiS  PAGE 


ttJjy  3C9_ 


1,  REPORT  security  classification 

unclassified 


REPORT  DOCUMENTATION  PAGE 


1t>  RESTRICTIVE  MARKINGS 


OlSiTV  CUSSiFiCaTion  AUTHORITY 


:d  OS  CLASSIFICATION -DOWNGRADING  SCMEOULE 


«  PERFORMING  ORGANIZATION  REFORT  NUMBERiS) 


3  OISTRI8UTION/AVAI  LABILITY  OF  REPORT 

Approved  for  public  release;  distribution 
unlimited 


5  MONITORING  ORGANIZATION  REPORT  NUM8E R(S) 

AFGL-TR-87-0229 


6«  NAVE  CF  PERFORMING  ORGANIZATION  |6b  OP  f  ‘CE  S  Y  MBO  L  7*.  NAME  OF  MONITORING  ORGANIZATION 

I  ilf  applicable ) 

Air  Force  Geophysics  Laboratory 


Saint  Louis  University 


oc  AO  ORE  S3  •  City .  Hate  and  /.IP  Cod*  » 

221  North  Grand  Blvd. 
St.  Louis,  MO  63103 


7b.  AOORESS  tCily.  State  and  ZIP  Coat) 

Hanscon  Air  Force  Base,  MA  01731 


Bb.  OFFICE  SYMBOL  9  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 
lit  applicable. 


Sc  ADDRESS  C.ty .  a  tale  and  ZIP  Code ) 


F19628-85-K-0029 


10  SOURCE  OF  FUNDING  NOS 


1400  Wilson  Blvd. 
Arlington,  VA  22209 


TASK 

WORK  UNiT 

NO 

NO 

DA 

AO 

13.  PERSONAL  AtTMORlSl 

R.  B.  Herrmann 


V  13*  T1?E  Cf  REPORT  13t.  TIME  COVERED  14  DATE  OF  REPORT  lYr  ,  Mo..  Day) 

l  Final  Report  from  1  APR  83 TO  31  MAR  87  87JUL15 


15.  PAGE  COUNT 

82 


18  SUBJECT  TERMS  (Continue  on  reverie  if  ncctuary  and  identify  by  block  number > 


Lg  wave  scattering,'  Phase  matched  surface-wave  filters 


'9  ABSTRACT  continue  on  reverie  if  neceiaorf  and  identify  by  bloc *  number 


fT'ne  scattering  of  Lg  waves  by  point  heterogeneities  in  a  crustal  waveguide  is 
considered.  Computer  simulations  lead  to  inferences  of  the  behavior  of  the 
(  scattering  mechanism  which  helps  understand  observed  facts  about  the  shear 
£  wave  coda,  especially  the  relationship  between  the  lag  time  at  which  measurements 
j'  are  made  and  the  inferred  coda  Q. 


A  second  study  reports  on  problems  in  some  routine  procei sing  techniques  used  for 
long  period  surface  wave  analysis.  The  influence  of  spectral , shape  on  inferred 


I  spectral  amplitudes  is  discussed. 


:0  OisTRiBuTiCN  AVAILABILITY  of  abstract 
4  UNCLASSIFIED  unlimited  □  SAME  AS  RPT.  □  OTIC  USERS  □ 


23«  NAME  OF  RESPONSIBLE  INDIVIDUAL 

James  F.  Lewkowicz 


DO  FORM  1473.  83  APR  edition  of  i  jan  73  >s  obsolete 


31  ABSTRACT  SECURITY  CLASSIFICATION 

Unclassified 


32b  TELEPHONE  NUMBER 
t/nclude  Area  Code t 

(617)  377-3028 


22l  OFFICE  SYMBOL 

AFGL/LWH 


SECURITY  CLASSIFICATION  OF  THiS  PAGE 


SECURITY  CLASSIFICATION  Of  Tmi*  PAGE 


lTlWlTO.  fMi A  Afl  V 


s."  R21  Mt  *«_•>  ha wtAiiA/< 


Table  of  Contents 


INTRODUCTION  1 

Synthesis  of  coda  waves  in  layered  medium  2 

Application  of  frequency  variable  filters  to  surface-wave  amplitude  analysis  41 


yyyyyvw w.%  ww  v  .v.v  v/.  uuvwt  v-  v  <jou  o  otv 


FINAL  TECHNICAL  REPORT 


Lg  Wave  Excitation  and  Propagation  in  the  Presence 
of  One-,  Two-  and  Three-Dimensional  Heterogeneities 


INTRODUCTION 

This  is  the  final  technical  report  for  this  contract.  Work  has  been  directed  both  toward 
understanding  seismic  wave  propagation  in  the  presence  of  simple  point  heterogeneities  in  the 
medium  and  also  toward  the  important  validation  of  computational  techniques  used  in 
analysis  and  synthesis  of  time  histories.  The  research  performed  has  been  submitted  for  pub¬ 
lication.  Two  papers  have  already  been  published.  The  papers  in  this  report  have  been  sub¬ 
mitted  for  publication.  The  paper  by  Wang  and  Herrmann  has  been  submitted  to  PAGEOPH 
for  a  special  volume  on  scattering  and  Q  edited  by  Aki  and  Wu.  The  paper  by  Russell, 
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Synthesis  of  Coda  Waves  in  Layered  Medium 

1  2 
Chien-Ying  Wang,  and  Robert  B.  Herrmann” 

^National  Central  University,  Department  of  Geophysics,  Chung-Li,  Taiwan,  Taiwan 

^Department  of  Earth  and  Atmospheric.  Sciences,  Saint  Louis  University,  PO  Box  8099,  St. 
Louis,  Missouri  63156,  USA 

Abstract  -  Malin’s  (1980)  first-order  single  scattering  theory  has  been  extended  to  study 
the  scattering  of  surface  waves  as  well  as  body  waves  by  distributed  point  scatterers  in  a  lay¬ 
ered  medium.  The  scattered  waveform  itself  is  generated  and  examined  instead  of  its  energy 
envelope.  The  theory  used  allows  1)  mode  conversion  2)  wave  type  conversion  3)  finite 
scatterer  distribution,  and  4)  the  effect  of  attenuation  from  scattering  as  well  as  intrinsic 
absorption.  The  cases  studied  are  for  elastic  or  slightly  attenuative  media  with  any  kind  of 
source  and  receiver  at  any  place  in  the  layered  structure.  This  direct  calculation  of  coda 
waves  provides  us  an  immediate  description  of  the  relation  of  coda  and  scattering.  The  objec¬ 
tives  are  to  find  1)  the  effect  of  layering  on  scattering,  2)  the  effect  of  scatterer  distribution  on 
recorded  vertical  and  horizontal  motion,  3)  the  relation  of  scattering  O  to  intrinsic  Q,  4)  the 
scattering  behavior  of  surface  and  body  waves,  and  5)  the  superposition  of  scattering  waves  to 
form  the  coda.  The  generation  of  body  waves  by  ’locked  mode’  approximation,  which  makes 
the  body-wave  scattering  a  subset  of  the  ’surface-wave’,  scattering  Preliminary  results 
explain  some  observed  coda  behavior  surprisingly  well.  NVe  find  a  larger  geometrical  spreading 
for  near  scatterers,  which  is  caused  by  mode  conversion  or  wavetype  conversion  because  of  the 
wide  angle  scattering.  This  makes  the  spreading  correction  higher  for  early  part  of  coda  which 
may  account  for  the  low  Q  observed  in  early  coda  of  regional  earthquakes.  This  study  is  of 
practical  value  as  an  effort  to  understand  the  complicated  coda  phases. 

Key  words:  Surface  waves,  scattering,  layered  media 

1.  Introduction 

Because  of  recent  interest  in  the  extraction  of  useful  information  from  the  coda,  a 
theoretical  study  of  coda  waves  is  urgently  required,  especially  since  the  theory  lags  behind 
observations.  Most  of  current  coda  analyzes  are  based  on  Aki’s  (1969)  backscattering  model. 
This  model  treats  the  coda  as  a  smoothly  decaying  envelope  which  is  formed  by  many  ran¬ 
domly  scattered  waves  Since  the  scattered  wave  is  the  result  of  an  averaging  process  for 
waves  reacting  with  randomly  distributed  scatterers,  most  theories  take  a  stochastic  approach 
and  describe  the  total  wavefield  by  several  statistical  quantities.  Mean  power  spectral  density 
is  frequently  used,  for  instance  (Aki,  1969;  Sato,  1977).  Using  such  an  approach,  the  analysis 
of  data  provides  a  gross  description  of  the  distributed  inhomogeneities  and  the  scattering  pro¬ 
cess,  but  some  ambiguities  arise.  The  problem  of  distinguishing  between  scattering  Q  and 
intrinsic  Q  is  a  typical  example  (.Aki,  1982,  Dainty,  1984).  For  our  study,  we  choose  a  deter¬ 
ministic  approach  to  investigate  the  effect  of  scattering  on  the  generation  of  coda.  We  will 
directly  examine  each  individual  scattered  pulse  generated  in  the  layered  medium. 

Sato  (1984)  proposed  a  scattering  model  which  describes  on  the  scatterers.  A  similar 
model  is  also  found  in  Wu  and  Aki  (1985).  This  model  considers  wave  type  conversion, 
scattering  pattern,  and  scatterer  size,  but  only  for  body  waves  in  a  homogeneous  w  hole  space. 
From  another  view  point,  Malin  (1980)  constructed  a  scattering  model  using  surface  wave  nor¬ 
mal  mode  theory.  An  important  feature  of  this  model  is  that  it  is  able  to  handle  multi-layered 
medium.  The  rays  scattered  by  the  inhomogeneitics  in  a  layered  structure,  such  as  the  earth's 
crust,  are  numerous.  It  is  only  by  using  mode  theory  that  scattering  in  such  a  medium  can  be 
treated . 

In  this  study,  we  extend  Malin’s  first-order  single  scattering  model  to  create  the  scatter¬ 
ing  ’pulse’  from  scatterers  situated  in  the  layered  structure.  A  simple  review  of  this  method  is 
provided  with  the  derivation  largely  simplified.  Each  individual  scattered  wave  is  examined 
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in  detail  rather  than  its  statistical  average.  The  study  is  directed  toward  understanding  the 
nature  of  the  scattered  waves  under  different  conditions.  This  enables  us  to  analyze  the 
behavior  of  each  scatterer,  and  hence  makes  it  easier  to  isolate  important  factors  which  may 
be  lost  in  an  averaging  process. 


2.  First-Order  Single  Scattering  Theory 

The  model  considered  is  a  layered  medium  (Fig.  1)  with  inhomogeneities  distributed  at 
some  restricted  regions  in  it.  The  inhomogeneous  region  is  assumed  to  be  small  compared  to 
the  wavelength.  We  call  it  a  scatterer  parcel.  If  the  parcel  is  not  small  enough  to  be  thought 
of  as  a  point,  we  can  divide  it  into  several  smaller  parts  and  perform  an  integration  to  find 
the  total  field.  For  this  study  we  just  consider  the  point  scatterer.  This  point  scatterer  gen¬ 
erates  secondary  waves  when  acted  upon  by  the  incident  wave.  The  scattered  wave  has  the 
form  of  pulse  with  limited  duration  (Fig.  2).  The  coda  is  thought  to  be  composed  of 
numerous  such  scattering  pulses. 

To  satisfy  the  Born  approximation  (weak,  single  scattering),  the  inhomogeneities  are 
assumed  to  be  small  perturbations  of  the  background  material: 

P  —  Po  +  fy 

Qjki  =  Co  ijki  +  K3ijkl 

p  is  the  density  and  Cyk|  is  the  elastic  constant.  For  the  isotropic  elastic  medium  considered 
here,  the  Cyk|’s  reduce  to  the  Lame’s  constant  X  or  the  rigidity  p.  Since  we  just  discuss  the 
scattering  from  a  small  parcel,  the  mean  behavior  of  these  inhomogeneities  is  not  assumed. 
These  three  parameters  p,  X,  and  p  can  be  reduced  to  a  single  parameter,  if  there  exists  a 
relation  among  them.  Sato  (1984)  used  the  velocity  perturbation: 

SVP  _  SVs 

vT"  =  e 

and  adopted  Birch’s  law  relating  the  density  and  wave  velocity 

Vp  =  ai  p  +  ao  Vs  =  b,  p  +  b0 


(a’s  and  b’s  are  constant).  This  results  in  , 

,  5X  bp  .  , 

-  =  vf  -T-  =  —  =  2  +  v)  ? 

p  X  p 

The  value  of  v  is  kept  at  0.8  in  this  report.  {  will  be  the  only  parameter  needed  to  describe 
the  inhomogeneity  variation  under  these  assumptions.  Wu  (1984)  used  different  values  for 
density  perturbation  and  elastic  constant  perturbation  which  results  in  two  different  types  of 
inhomogeneities;  impedance  type  and  velocity  type,  each  of  which  has  different  scattering  pro¬ 
perties.  To  reduce  the  number  of  independent  parameters,  we  just  use  one  parameter  (  to 
represent  them. 

As  indicated  by  Hudson  (1977),  the  scattered  field  for  any  first-order  scattering  theory 
can  always  be  expressed  using  the  representation  theorem: 


uj"  =  JdV 


bp  < J 2  GiJ  Ui°  -  «5Cyk,  dp,1  a,uk° 


For  the  isotropic  elastic  case,  this  is  equivalent  to 

uj')  =  JdV  j  bp  u 2  G;J  Ui°  -  b\  dp1  djUj°  -  bp  dfi1  (dp?  +  d,u°)  J 


(1) 


-  4  - 


I 


where  u°  is  the  wave  field  incident  on  the  inhomogeneous  body,  and  GjJ  is  the  Green’s  func¬ 
tion  for  the  corresponding  boundary  conditions.  For  the  case  of  body  wave  in  a  homogeneous 
space,  GiJ  has  the  form  as  given  by  equation  (4.30)  in  Aki  and  Richards  (1980),  which  is  used 
by  Sato  (1984)  and  Wu  and  Aki  (1985).  For  the  case  of  a  surface  wave  in  the  layered 
medium,  GjJ  can  be  determined  by  the  normal  mode  theory  (Malin,  1978). 

The  surface  wave  case  is  easy  to  construct  since  the  incident  wave  and  the  scattered 
wave  are  trapped  in  the  layers  and  propagate  two  dimensionally.  From  surface  wave  mode 
theory  (Haskell,  1963;  Haskell,  1964;  Wang,  1981),  the  incident  wave  field,  which  has  mode 
order  m  and  wave  type  v  (  v  =  R  for  Rayleigh  wave  and  L  for  Love  wave),  has  the  form 


"am  = 


V27T  l'kmx 


-i 

"Sn,  e  4  ‘Vn 


o  V  v  X 

m  om  *0m 


*U  =  [  Ur  ,  0  .  -iU, 
LU  =  [  0  ,  U,  ,  0 


at  scatterer 


(2) 


where 

k  =  wave  number, 

x  =  the  distance  between  the  source  and  the  scatterer, 
a  =  amplitude  factor, 

C  =  phase  velocity, 
g  =  group  velocity, 

Io  =  energy  integral  of  surface  waves. 

Ur,  U„  U|  =  displacement  eigenfunctions  in  the  radial, 
vertical  and  tangential  direction, 


The  source  term  Sm  for  double-couple  and  explosive  sources  is  given  in  the  Appendix.  Solv¬ 
ing  equation  (1),  the  Green’s  function  with  mode  n  and  wave  type  rj  for  the  wave  field  after 
scattering  is 

i  -i  "k-r-t-i- 

"GnJ  =  —J- -  TJd  X  e  4  ”Rdj  (3) 

V2jt  "k.r 


Ru  = 

y  Urcos0  ,  U,sin<i  ,  —  iU, 

J  at  scatterer 

"R  = 

Urcos0  ,  Ursin<^  ,  —  iU, 

.  at  receiver 

LU  = 

Ufsinff  ,  —  U<cos0 , 0  j 

at  scatterer 

LR  = 

y  U#sin<l>  ,  — U<cos <t>  ,  0  j 

at  receiver 

r  is  the  distance  from  the  scatterer  to  the  receiver.  The  scattering  angle  9  and  the  receiving 
angle  <t>  are  defined  in  Fig.  1.  These  angles  are  important  since  they  control  the  recorded 
wavefield.  Scattered  P-SV  waves  can  appear  on  the  receiver  horizontal  component  defined  as 
tangential  with  respect  to  the  source 

Substituting  these  functions  into  the  equation  (1),  the  scattered  wave  field  which 

includes  the  conversion  from  mode  m  of  wave  type  v  to  mode  n  of  wave  type  v.  is  written  as 


"'ul’i-ZdV 
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2s-V  *Tcmx  '7k„r 


am  ”ans, 


V)+ 
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(tp“2  "’A ma-S\  ”KmD-6n  ”>Cma)  (4) 

where 

Amo  =  ^rm  U  a  cos#  +  Ulm  Lln 
RRBm  =  (  Rkm  Rurm  -  d,u,m)  (  Rka  Urn  -  32U2a) 

""C™.  =  2  (  BkmUrm  Rk„Urncos2#  +  d2U2(nd2U2a) 

+  ( Rkmu,m  +  a2urm)  ( Rkau2n  +  a2urn)  cos# 

^Amo  =  -Urm  U,B  sin# 

“Cm.  =  ~  RkmUrm  LkaU,a  sin  2# 

-a2U,0  (3,Urm  +  RkmUjm)  sin# 

^A™  =  U,m  Urn  sin# 

HcmU#m  RkaUrn  sin2# 

+  dtV,m  (d2Ura  +  RknU2n)  sin# 


^Amj  =  Ufa,  Ufa  COS# 

UCmn  =  ^mUlm  ^Ufa  COS2#  +  <9jU|m  d,Ufa  COS# 


We  also  define  a  scattering  term  Fmn  as 
-i(  *1tmx+ "k.rj+if 


^m.  -  fdV 
=  JdV 
=  JdV 


(W*  ^Amn-#X  "’Bron-#/t  "'CBJ 


-Orf+V>*f 


«  [v/w2  ^AmD-(2+v)X  "'Bm,-(2+v)/i  "’Cmn] 


-i(  »krf+V)-Hv 


(5) 


The  final  total  field  solution  is  the  sum  over  all  modes  and  conversions  for  all  scatterers  in  the 
region.  The  scattering  term  F  corresponds  to  the  ’volume  factor’  of  Wu  and  Aki  (1985)  if 
the  Green’s  function  related  term  D  is  not  included.  For  a  point  scatterer,  we  simply  use 


"Tm.  =  SV  e 


-i(  l’kmx+  ’fknr)+i’ 


?  "?D„ 


with  SV  having  the  dimension  of  unit  volume. 


S.  Scattering  Q 

To  include  the  effect  of  energy  loss  due  to  scattering,  Malin  (1978)  considered  a  correc¬ 
tion  which  conserves  energy  to  the  first  order  by  choosing  a  specific  correlation  function  for 
the  scattering  cross  section.  Although  we  are  interested  in  the  behavior  of  each  individual 
point  scatterer,  the  scattering  coefficient,  or  equivalently  the  cross  section,  is  still  needed  when 
defining  the  scattering  Q. 

Let  us  consider  a  small  scatterer  parcel  containing  distributed  inhomogeneities.  For  a 


single  point  scatterer,  we  can  use  a  deterministic  scattering  formula.  However,  when  a  wave 
is  incident  on  this  scatterer  parcel,  the  size  of  the  scatterers  within  it,  the  relative  distance 
between  individual  point  scatterers,  and  the  wavelength  of  the  incident  wave  will  determine 
the  partition  of  energy  being  scattered.  From  scattering  theory,  the  parameter  describing  this 
property  is  the  ’cross  section’  (Ishimaru,  1977,  p.10  ).  Within  a  unit  of  the  scattering  angle, 
the  differential  cross  section  crd  is  defined  by 


..  r  (  scattered  energv  flux  density) 

~  llm  — ,■ — ^ : - : — t-2 

r— oo  (  incident  energy  flux  density) 


The  total  scattering  cross  section  is  obtained  by  integration  over  all  angles: 

<r,  =  ferd  d& 

2* 

When  scatterers  are  distributed  with  the  mean  density  n  (particles  per  unit  area  at  a  constant 
depth  for  a  two  dimensional  case),  the  scattering  coefficient  a  which  characterizes  the  inten¬ 
sity  of  the  scattered  wave  excitation  is  given  by 

a  =  n  cr, 


The  reciprocal  of  a  is  usually  called  the  mean  free  path.  The  scattering  coefficient  a  thus 
defined  is  equivalent  to  the  turbidity  in  the  two  dimensional  case  when  compared  to  Aki  and 
Chouet  (1975). 

To  obtain  the  energy  flux  density,  the  wave  field  excited  in  the  layered  medium  should 
be  taken  into  account: 

<  | u  |2>  =  w2  Jdzpo  < uu* > 


<>  means  the  average  over  the  ensemble.  Using  equations  (2)  and  (4),  it  is  not  difficult  to 
find  the  scattering  coefficient  for  mode  pairs  mn,  and  types  urj  by  setting  r  in  equation  (6)  as 
the  distance  between  the  scatterer  and  the  receiver: 


">a 


mo 


=  n  Jdff 
2* 


r  <  I  ^)|2> 

<  I  "u°|2> 


2tt 


W  %■> 

nK  "Ion, 


/ "  <  "T2 


>  d0 


/ n  <  "T2^  d 0 


/ n  <  "'F2n>  d9 

2n 


Now  we  have  the  problem  of  defining  the  average  behavior  of  the  function  <  ‘/,?F2n>  . 
Taking  a  closer  look  at  the  scattering  region,  Fig.  1,  the  scatterers  are  distributed  over  a  plane 
area  and  the  phase  term  in  F£n  (equation  5)  is  actually 

_  e-'«'cdol<kn?‘x-knrr) 

_  e_ ‘km*~ 'kor  e'T  k*B 

ke  =  k„— km  is  called  the  Bragg  vector.  Hence,  in  evaluating  FroD  from  equation  (5), 
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where  we  have  made  a  two  dimensional  Fourier  transform  for  the  scatterer  distribution  spec¬ 
tra  in  the  horizontal  direction  (Kennett,  1072).  The  ensemble  average  of  Fmn  is  thus 


<F2n>  =  Jdzdz 


<  C(kB.z)  )  >  Dm0(z)  ) 


<££  >  ,  the  correlation  function  of  inhomogeneity  distribution,  is  usually  given  by  a  form  of 
an  exponential  function  or  a  Gaussian  function  (Chernov,  1960).  Assuming  that  this  function 
has  a  small  side  lobe,  i.e.,  small  correlation  distance  in  the  exponential  form,  we  have  a 
simpler  form: 

n<F2n>  =  [  n Jdzdz'  < £(kB,z)  £(kB,z  )>  ]  D*n(z)  =  D2a(z) 


The  term  in  the  bracket,  which  describes  the  condition  of  scatterer  distribution  around  the 
point  considered,  will  be  kept  together  and  called  a  scattering  attenuation  unit  (Usa)  .  U„  has 
the  units  of  km4.  If  USI  is  maintained  as  a  constant  value  in  the  calculations,  this  is 
equivalent  to  stating  that  the  scattering  environment  is  similar  for  all  waves  propagating  in 
the  layered  medium.  The  expression  above  is  perhaps  too  simplified.  The  correlation  func¬ 
tion  must  be  specified  if  an  accurate  scattering  cross  section  is  desired.  In  this  study,  we  just 
use  the  simplest  form  to  examine  its  possible  effect. 


The  scattering  coefficient  thus  has  the  from: 
""a  =  mAU„  +  I"K©2+  «*e0 
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0’s  are  given  in  the  Appendix  I.  We  have  used  the  distance  r  between  the  scatterer  and  the 
receiver  to  define  the  scattering  coefficient  (equation  6)  and  to  calculate  the  incident  energy 
flux  and  the  scattered  energy  flux  by  fixing  the  source-scatterer-receiver  geometry.  The  mode 
pair  mn  reduces  its  energy  by  e  mn  after  traveling  the  distance  from  source  to  scatterer  to 
receiver.  The  effect  thus  caused  by  the  scattering  of  inhomgeneities  is  called  the  scattering 
attenuation. 


4-  Numerical  Experiment  I:  Wave  type  conversion  and  Mode  conversion 

With  the  velocity  model  chosen  (Cl’S  model,  see  Table  1),  we  generate  the  eigenfunc¬ 
tions  of  surface  waves,  mainly  the  Lg  wave,  up  to  5  liz.  Using  these,  the  scattered  waves  are 
created  by  specifying  the  source  (depth,  mechanism),  the  receiver  (distance,  instrument),  the 
scattcrers  (location,  material  variation,  and  scattering  intensity),  and  the  attenuation 


condition  (intrinsic  Q  and  scattering  Q).  To  make  the  problem  easier  to  follow,  we  fix  the  fol¬ 
lowing,  except  as  indicated  otherwise: 

source  depth  —  10  km 
scatterer  depth  =  -1  km 

source  time  function  =  parabolic  with  base  0.4  second, 
receiver  distance  —  100  km 
waveform  observed  =  ground  velocity 
inhomogeneity  variation  —  5% 

mode  conversion  =  up  to  10  neighboring  modes  (due  to  computer 
restriction  on  speed  and  storage) 
wave  type  conversion  =-  yes 

An  arbitrary  source  mechanism  was  used  and  the  effect  of  its  radiation  pattern  is  averaged 
out  by  distributing  several  scatterers  on  the  ellipse.  In  this  section,  we  will  discuss  the  proper¬ 
ties  of  t  e  scattered  waves  under  different  scattering  conditions  but  ignore  the  problem  of 
scattering  and  intrinsic  attenuation  at  this  moment.  The  so  called  ’diagonal  selection’  rule 
(Malin,  1980)  which  states  that  a  particular  mode  of  surface  wave  mainly  scatters  into  the 
same  mode  without  significant  conversion  will  be  examined  first.  Scattering  from  different 
locations,  e.g.,  distances,  depths,  scattering  angle,  etc.,  are  also  discussed. 

Malin  (1980)  used  the  scattering  coefficient  of  vertical  component  from  an  acoustic 
model  to  justify  the  diagonal  selection  rule.  Since  the  scattered  wave  instead  of  its  energy 
envelope  is  generated  here,  this  rule  will  be  discussed  by  examining  the  resultant  waveforms. 
Two  cases  are  studied;  one  is  for  the  wave  type  conversion  (from  R  to  L  and  L  to  R),  and  the 
other  for  the  mode  conversion.  Fig.  2  shows  the  waveforms  of  signals  directly  from  the  source 
(  surface  wave)  or  from  a  scatterer  (scattered  waves)  after  traveling  the  same  distances  as 
indicated  at  the  ends  of  seismograms.  Both  the  source  depth  and  the  scatterer  depth  are  kept 
at  10  km.  Scattered  waves  have  more  high  frequencies  and  are  of  shorter  duration  than  the 
direct  surface  waves.  The  major  difference  is  in  the  excitation  of  the  fundamental  mode. 
The  fundamental  mode  excited  by  the  source  does  not  excite  the  scatterers  as  well  at  high  fre¬ 
quencies  as  do  the  higher  modes.  This  is  a  result  of  the  eigenfunction  distribution  with  depth. 

Fig.  3  examines  the  effect  of  wave  type  conversion  on  scattering.  In  this  figure  and 
those  that  follow,  the  number  of  kilometers  indicates  the  total  scattering  distance,  i.e.  the  dis¬ 
tance  from  the  source  to  the  scatterer  plus  the  distance  from  the  scatterer  to  the  receiver.  The 
distance  from  the  source  to  the  receiver  is  fixed  at  100  km.  The  angle  in  degrees  is  the  loca¬ 
tion  of  the  scatterer  on  an  ellipse  which  has  the  source  and  the  receiver  at  its  foci  and  the 
coordinate  origin  at  its  center.  Zero  degrees  indicates  that  the  scatterer  is  behind  the  receiver, 
and  90  degrees  indicate  that  the  scatterer  is  equidistant  from  the  source  and  receiver.  Three 
different  components  of  motion  are  plotted,  vertical  Z,  radial  R  and  the  tangential  T  with 
respect  to  the  source-receiver  coordinates.  In  each  pair  of  traces  of  Fig.  3,  the  lower  one 
includes  the  wave  type  conversion  and  the  upper  one  does  not.  No  mode  conversion  is 
involved.  From  this  comparison,  we  find  that  the  wave  type  conversion  is  not  obvious  except 
at  short  distances  for  wide  angle  scattering.  The  Z  component  seems  to  be  more  independent 
of  this  effect  than  the  other  components.  Thus  it  is  concluded  that  wave  type  conversion  is 
important  only  for  the  horizontal  components  at  short  distances.  It  is  at  these  distances  that 
wide  angle  scattering  occurs.  For  large  scattering  distances  relative  to  the  source-receiver  dis¬ 
tance,  backscattering  is  all  that  occurs  from  the  geometry. 

Fig.  4  illustrates  the  effect  of  mode  conversion.  Three  cases  are  displayed  in  which  the 
number  of  modes  permitted  to  be  converted  to  other  near  neighboring  modes  are  1,  5  and  10 
In  addition,  the  effects  of  wavetype  conversion  are  also  shown.  The  scattered  waveforms  are 
created  by  combining  the  effect  of  seven  scatterers  evenly  distributed  on  the  ellipse  defined  by 
the  scattering  distances.  The  amplitude  scale  of  plotting  for  every  scattering  distance  is  the 
same.  As  in  Fig.  3,  the  Z  component  seems  stable  and  does  not  accumulate  as  much  contribu¬ 
tion  from  other  modes  except  at  the  short  distance  where  mode  conversion  must  still  be 


considered.  On  the  other  hand,  the  R  and  T  components,  Fig.  4(b)  and  4(c),  are  sensitive  to 
the  number  of  modes  allowed  to  convert  each  other  and  to  the  wave  type  conversion.  These 
results  do  not  support  the  diagonal  selection  rule,  although  the  use  of  it  alone  will  not  cause 
major  differences  in  the  coda  excitation.  For  an  elastic  medium,  the  energy  exchange  by  mode 
conversion  is  important  especially  for  the  two  horizontal  components  and  for  short  distances. 
The  mode  conversion  can  be  looked  upon  as  being  similar  to  a  change  of  incident  angle  on 
the  layered  medium  after  scattering.  This  is  an  effect  similar  to  the  scattering  pattern  for 
body  wave  illustrated  in  Aki  and  Richards  (1980,  Chap  13).  At  short  scatterer  distances, 
because  of  the  variation  of  scattered  angle,  the  mode  conversion  must  play  an  important  role, 
affecting  the  energy  of  scattered  waves. 

We  next  examine  the  effect  at  different  scatterer  depths.  A  result  is  shown  in  Fig.  5  for  a 
scattering  distance  of  200  km  and  a  source  depth  of  25  km..  A  Q  model  10U0.5Q,  which  will 
be  defined  later,  is  also  applied.  It  is  found  that  the  scatterer  at  greater  depth  preferentially 
excites  the  scattered  waves  with  the  higher  phase  velocities.  The  scatterer  in  the  upper  sedi¬ 
mentary  layer  (0.5  km)  is  capable  of  generating  well  dispersed,  long  duration  scattered  waves. 
Because  of  the  dissipative  attenuation  included,  the  shallow  scatterers  yield  waveforms  that 
lose  their  energy  and  frequency  content  faster  with  increasing  distance,  which  can  make  them 
similar  to  the  waveforms  of  deeper  scatterers.  At  larger  propagation  distances,  the  difference 
between  scattered  waves  at  different  depths  is  gradually  reduced,  although  the  effect  of 
attenuation  may  still  be  seen.  Our  studies  indicate  that  a  scatterer  at  any  depth  in  the  crust 
generates  a  similar  scattered  signal,  as  long  as  the  propagation  distance  in  large. 

5.  Numerical  experiment  II:  Scattering  Q  and  Intrinsic  Q 

As  shown  in  the  previous  section,  the  scattering  coefficient  is  a  function  of  the  correla¬ 
tion  between  scatterers,  the  inhomogeneity  size,  and  the  incident  wave  length.  To  make  the 
problem  easier  to  solve,  we  defined  a  scattering  attenuation  unit  (U,J  from  which  we  obtained 
the  scattering  coefficient  a  after  taking  into  account  the  angular  scattering  variation.  The 
amplitude  decay  is  given  by  e-0^2  where  r  is  the  total  scattering  distance.  In  fact,  this  is  a 
definition  of  turbidity  which  is  used  to  describe  the  scattering  energy  loss  mechanism  (Dainty, 
1981).  In  our  model,  the  scattering  attenuation  is  assumed  the  same  for  different  scatterers 
which  are  distributed  randomly  and  uniformly  along  scattering  wave  path  passing  through 
them.  We  will  examine  the  effect  of  scattering  on  the  wave  energy  loss  and  compare  this  with 
the  attenuation  due  to  anelastic  absorption.  The  anelastic  attenuation  is  calculated  by  per¬ 
turbing  the  elastic  constants  and  combining  with  the  intrinsic  Q  values  using  surface-wave 
variational  techniques. 

Fig.  6  shows  an  example  with  scattering  distance  at  300  km  and  a  scatterer  depth  of 
0.5km.  The  uppermost  trace  in  this  figure  is  the  reference  trace  which  does  not  involve  any  Q 
effect.  Other  traces  are  named  by  the  symbol  pl'qQ  which  is  read  as  scattering  Q  having  p 
U„  and  intrinsic  Q  having  values  from  the  model  q  The  basic  Q  model  used  is  Qy  -■  300  for 
upper  20  km  and  2000  for  the  rest  in  the  Cl’S  model.  20  indicates  twice  as  much  attenua¬ 
tion  as  IQ,  eg.,  Q^=150.  It  is  obvious  that  the  scattering  Q  reduces  the  high  frequency 
energy  less  than  intrinsic  Q.  This  means  that  the  scattering  Q  is  not  too  sensitive  to  the 
scattering  environment.  If  the  anelastic  part  of  structure  does  not  dissipate  all  high  frequency 
signal,  the  scattering  Q  may  not  reduce  it  even  under  strong  scattering  conditions  Further¬ 
more  the  scattering  Q  reduces  the  scattered  pulse  amplitude  but  does  not  affect  its  frequency 
content  much  as  can  be  seen  from  the  traces  in  101  0Q  and  lOl’OO.  We  also  find  that  the 
scattering  Q  seems  to  suppress  the  lower  mode  signals  and  uihanee  the  high  frequency  higher 
modes.  However,  this  feature  probably  depends  on  the  velocity  model  used  From  this  test, 
we  conclude  that  the  effect  of  scattering  Q  and  intrinsic  Q  are  apparently  different 

If  the  scattered  arrivals  from  different  distances  constitute  the  coda  at  different  lapse 
times,  the  frequency  content  of  these  pulses  may  indicate  the  effect  of  O  at  different  places  of 
coda  To  examine  this  effect,  we  generate  the  scattered  pulses  at  a  distance  of  500  km  under 


different  Q  conditions  These  are  shown  in  Fig  7  The  Parts  a  and  b  are  for  the  scatterer 
depths  at  5  km  and  15  km  at  these  depths  the  values  are  different.  The  spectra  of  the 
traces  are  also  plotted,  shifted  by  a  factor  of  ten  for  clarity.  One  thing  revealed  from  this  test 
at  different  distances  is  that  the  scattering  attenuation  absorbs  the  high  frequency  signal 
equally  well  irrespective  of  the  scattering  distance  or  the  scatterer  depth.  This  is  due  to  the 
fact  that  a  similar  scattering  environment  is  assumed  when  calculating  scattering  coefficient  a 
The  intrinsic  Q,  on  the  other  hand,  attenuates  the  high  frequency  signals  rapidly  for  a  shal¬ 
lower  scatterer  where  intrinsic  Q  is  low.  Fig  7(a)  and  7(b)  can  be  thought  of  as  two  extreme 
cases:  in  7(a)  intrinsic  Q  dominates,  and  in  7(b)  scattering  Q  dominates.  The  final  Q  value 
must  result  from  the  interaction  of  these  two  attenuation  mechanisms.  If  the  data  show  a 
strong  frequency  dependence  at  high  frequency  range,  it  is  possible  this  is  the  attenuation 
from  deep  structure  and  the  effect  of  scattering  Q  prevails.  At  the  low  frequency  (<1  Hz),  it 
is  difficult  to  tell  the  difference  between  these  two  mechanisms. 

6.  Numerical  Experiment  III:  Geometrical  Spreading 

Aki  (1969)  proposed  a  scattering  theory  which  describes  the  coda  waves  of  a  local  event 
as  being  composed  of  backscattered  waves  from  many  random  distributed  inhomogeneities  in 
the  lithosphere.  A  relation  from  his  model,  which  has  been  widely  used  in  the  literature,  is 

A(w  |t)  =  C(ui)t-‘'e^',/2Qc 

where  A(w  |t)  is  the  average  peak-to-peak  amplitude  around  the  frequency  w  and  the  time  t, 
C(w)  is  the  coda  shape  function  which  includes  the  scattered  wave  excitation  term  and  a 
dispersion  correction,  t-1/  is  for  geometrical  spreading,  and  the  exponential  term  describes  the 
dissipative  attenuation.  The  value  of  it  is  taken  as  1.0  for  body  waves,  0.5  for  surface  waves, 
and  0.75  for  a  diffusion  model  (Aki  and  Chouet  1975).  Treating  coda  of  local  earthquakes 
(distance  <  100  km),  the  source  and  the  receiver  are  set  at  the  same  place  to  simplify  the  cal¬ 
culation  of  backscattering  time  in  this  model.  Sato  (1977)  proposed  a  modification  of  the 
geometrical  spreading  factor  for  near  scatterers  which  create  waves  arriving  at  the  receiver 
just  after  the  direct  wave.  This  modification  takes  into  account  the  separation  between  the 
source  and  the  receiver  and  is  used  for  regional  earthquakes  (Pulli,  1984).  Kopnichev  (1977) 
also  derived  a  similar  correction  by  considering  the  ellipticity  of  the  scatterer  surface  for 
arrivals  at  a  given  time.  In  fact,  this  correction  represents  a  geometrical  change  of  the 
scatterer  distribution  ranging  from  a  flattened  ellipsoid  to  a  shape  similar  to  a  sphere  at  large 
scatterer  distances  (from  an  ellipse  to  a  circle  for  surface  waves  as  illustrated  in  Fig.  1.  It  is 
purely  a  geometrical  adjustment.  However,  as  discussed  in  preceding  sections,  the  scatterers 
at  the  near  distances  yield  significant  mode  conversion  and  wave  type  conversion.  A  more 
complicated  nature  of  scattered  waves  is  expected  there.  Similarly,  because  of  the  importance 
of  scattering  angle  in  controlling  the  scattering  mechanism  (Wu,  1984),  its  effect  must  be 
significant  at  short  distances  where  the  scattering  angle  changes  abruptly.  In  this  section,  we 
will  measure  this  effect  from  simple  numerical  studies.  The  result  indicates  that  the  so  called 
’Sato’s  correction'  which  gives  a  heavier  geometrical  spreading  correction  in  the  early  part  of 
coda  than  Aki's  formula  is  still  too  small  as  the  scattered  wave  energy  changes  rapidly 
between  modes  and  wavetypes  at  near  distances. 

In  order  to  see  the  variation  of  geometrical  spreading  factor,  we  generate  many  scattered 
pulses  at  different  scattering  distances,  ignoring  dissipative  attenuation,  both  of  intrinsic  and 
scattering.  (In  fact  four  distributed  scatters  are  used  for  each  scattering  distance,  and  result¬ 
ing  amplitudes  are  increased  by  the  square  root  of  the  distance  to  approximate  a  uniform  dis¬ 
tribution  of  scatterers).  An  averaged  peak-to-peak  amplitude  is  then  measured,  with  narrow 
band  filtering  if  needed,  and  plotted  on  log-log  scale  with  the  distance.  One  of  the  results  is 
shown  in  Fig.  8(a).  The  scattering  distance  is  spaced  by  10  km  up  to  GOO  km  and  the  receiver 
is  kept  at  100  km  from  the  source.  It  is  obvious  that  there  exist  two  slopes.  A  steeper  slope 
with  amplitude  dropping  rapidly  with  the  distance  is  seen  in  the  front  part.  A  similar 
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tendency  is  also  observed  for  the  tangential  component,  Fig.  8(b).  In  these  figures,  the 
unfiltered  scattered  wave  (upper  one)  and  the  filtered  scattered  wave  all  follow  similar  trends. 
This  is  due  to  a  similar  frequency  response  for  scattered  wave  even  traveling  over  different  dis¬ 
tances. 

Fig.  9  shows  a  similar  amplitude  decay  for  a  western  United  States  model  (Table  2). 
Here  we  have  set  the  source-receiver  distance  at  300  km,  Fig.  9(a),  and  at  100  km,  Fig.  9(b). 
The  rapid  amplitude  decay  happens  at  the  scatterer  distances  smaller  than  about  600  km  and 
200  km  respectively.  Hence,  we  may  define  an  ’unstable  region’  corresponding  to  scattering 
distances  smaller  than  twice  the  source-receiver  distance,  which  has  a  rapid  decrease  in  ampli¬ 
tude,  and  a  ’stable  region’  at  larger  distances  which  has  a  spreading  factor  close  to  unity. 
This  distance  separation  agrees  with  that  of  ’Sato  correction’  of  twice  the  source-receiver  dis¬ 
tance.  Note  that  the  decay  exponent  in  unstable  region  is  greater  than  that  in  stable  region 
by  0.5- 1.0.  This  value  is  larger  in  WUS  model  than  in  CUS  model  (Table  1)  for  the  Z  com¬ 
ponent.  Taking  this  into  account,  we  propose  a  geometrical  spreading  correction  to  the  RMS 
coda  amplitude  based  on  an  extension  of  Sato’s  (1977)  model  to  the  surface  wave: 


|2*(— +1) 

o&  2‘ 


for  a  <  2 


for  a  >2 


where  a  is  the  ratio  of  lapse  time  on  the  coda  to  the  Lg  wave  propagation  time  and  a  is 
between  0.5  and  1.0.  At  larger  lapse  times,  the  geometrical  spreading  correction  is  propor¬ 
tional  to  a-1.  This  correction  is  plotted  in  Fig.  10.  The  factor  involving  the  square  roots  is 
obtained  by  an  extension  to  the  body-wave  derivation  of  Sato  (1977).  The  term  following 
represents  a  correction  based  on  the  numerical  experiments. 

In  propagating  in  a  layered  medium,  not  only  does  the  amplitude  of  the  scattered  pulses 
decay  with  the  distance  but  its  duration  also  increases  because  of  surface  wave  dispersion.  We 
generated  a  set  of  waveforms  at  scattering  distances  of  110  to  450  km,  measured  the  signal 
duration.  Fig.  11  is  a  linear  regression  which  shows  the  increase  of  measured  duration  with 
propagation  distance.  The  CUS  model  exhibits  an  approximate  increase  in  duration  of  0.05 
sec/km  while  WUS  model  exhibits  an  increase  of  0.5  sec/km.  Because  a  thicker  sedimentary 
layer,  the  WUS  model  is  to  be  more  dispersive.  This  dispersion  induced  pulse  duration  expan¬ 
sion  must  be  considered  when  summing  all  pulses  to  form  the  coda. 


7.  Numerical  Experiment  IV:  Superposition  of  Scattered  Pulses 

We  have  discussed  the  properties  of  a  single  scattered  pulse.  The  next  question  is  how 
these  pulses  are  combined  to  form  the  coda.  In  this  section,  we  will  conduct  a  numerical 
experiment  using  the  parameters  determined  in  the  preceding  sections  to  synthesize  the  coda 
envelope.  From  recent  studies  of  coda  Q  of  regional  earthquakes,  Pulli  (1984)  and  Shin  (1985) 
both  observed  lower  coda  Q  values  from  the  early  part  of  coda.  Der  et  al.  (1984)  also  men¬ 
tioned  an  obvious  coherence  difference  in  early  and  later  coda.  This  section  will  attempt  to 
propose  an  explanation  to  these  observations. 

From  the  study  of  scattered  pulses  just  discussed,  we  make  the  following  assumptions; 
1)  the  shape  of  scattered  pulse  can  be  approximated  by  an  exponential  function  te~il;  2)  the 
width  of  the  pulse  depends  on  its  traveling  distance  with  duration  =  a0*distance;  3)  the  max¬ 
imum  amplitude  of  the  pulse  decays  as  e  t“m,  where  the  exponential  term  describes  the 
effect  of  attenuation  including  intrinsic  or  scattering  Q,  and  t-m  is  the  geometrical  spreading 
factor;  4)  the  pulses  arrive  at  the  receiver  in  a  random  manner,  thus  the  number  of  pulses  in  a 
time  interval  at  any  time  of  seismogram  is  assumed  to  be  approximately  constant.  To  take 
into  account  reasonable  variation,  we  also  allow  a  10fc  random  variation  for  the  parameters: 
m  (spreading  factor),  a,,  (duration  expansion)  and  dtp  (pulse  interval).  We  will  discuss  the 
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pulse  at  1  Hz  with  the  propagation  velocity  3  5  km/sec. 

With  each  pulse  defined,  the  coda  is  then  synthesized  by  summation  of  many  pulses. 
One  example  of  synthetic  coda  thus  formed  can  be  seen  from  Fig  12  where  every  individual 
pulse  and  the  coda  envelope  after  summation  are  plotted  assuming  different  parameters. 
Since  the  energy  of  random  arrivals  is  mixed,  we  use  a  root-mean-square  rule  to  form  the  coda 
envelope.  After  constructing  the  coda,  we  then  apply  the  Aki-Chouet’s  (1975)  model  with 
surface-wave  spreading  correction,  i.e.,  t“°  5  to  calculate  the  coda  Q  value.  This  corresponds 
to  an  m  =  1.0  for  the  geometrical  spreading  factor  of  the  individual  scattering  pulses.  The 
objective  here  is  to  see  how  the  coda  Q  estimate  is  affected  if  the  actual  geometrical  spreading 
in  the  early  part  of  the  coda  is  actually  more  rapid.  The  Q  values  at  two  regions,  first  100 
seconds  and  the  rest,  are  calculated  separately.  Here  a  sliding  window  with  10  seconds  length 
(Pulli,  1984)  is  also  used  when  taking  points  from  the  coda  shape.  The  cases  for  different 
values  of  m,  ao  and  dtp,  which  represent  different  scattering  conditions,  are  examined. 

One  interesting  result  revealed  by  this  simple  test  is  that  the  pulse  spreading  factor,  m, 
is  the  most  important  factor  which  controls  the  coda  shape  decay.  The  duration  factor  and 
the  pulse  density,  on  the  other  hand,  only  have  limited  effect  on  the  coda  formation  if  the 
pulses  are  crowded  together  enough  to  make  a  smoothly  decaying  coda  shape  as  usually 
observed.  The  dispersion  of  scattered  pulses  in  the  layered  medium,  which  causes  the  dura¬ 
tion  increase,  does  not  affect  the  coda  as  much  as  first  expected  (Dainty,  1984).  We  also 
found  that  for  a  low  Q  region  the  coda  decays  rapidly,  thus  not  too  many  pulses  are  needed 
and  its  Q  value  comes  out  fairly  stably.  The  coda  is  said  to  be  ’saturated’  by  attenuation. 
For  a  high  Q  region,  different  pulse  parameter  will  play  an  important  role  on  the  coda  forma¬ 
tion. 

As  a  result  of  repeating  the  above  experiment  for  different  combination  of  parameters  m, 
a,,,  and  dtp,  we  obtained  a  preliminary  understanding  of  the  constitution  of  coda  by  the  sum¬ 
mation  of  scattering  waves: 

1)  For  a  high  Q  region,  such  as  the  eastern  United  States,  the  manner  the  scattered  pulse 
decays  with  the  distance,  because  of  geometric  spreading  or  mode  conversion,  controls  the 
value  of  Q  determined  from  coda  shape.  A  large  m  value  in  the  early  part  of  coda  (unstable 
region)  has  a  major  efTect  on  Q  value.  This  is  the  case  predicted  in  Fig.  12(a). 

2)  For  a  low  Q  region,  such  as  the  western  United  States,  the  pulses  which  arrive  at  the 
receiver  to  form  the  coda  are  strongly  affected  by  the  attenuation  and  thus  the  coda  shape  is  a 
stable  estimator  of  Q.  This  is  show  n  in  Fig.  12(b). 

3)  The  pulses  with  longer  duration  yield  a  longer  tail  which  in  turn  leads  to  a  higher  Q  esti¬ 
mate,  but  this  effect  is  not  as  strong  as  expected. 

4)  If  the  pulses  do  mix  to  form  the  coda,  the  density  of  pulse  arrivals  is  not  very  important 
with  respect  to  the  final  coda  Q  values.  This  is  due  to  the  fact  that  the  pulses  are  superim¬ 
posed  by  means  of  the  energy  which  reduces  the  difference  in  the  number  of  pulses. 

5)  Estimation  of  coda  Q  using  only  earlier  parts  of  the  coda  requires  knowledge  of  the  rate 
of  correct  geometrical  spreading  term  there 

The  summation  of  pulses  as  mentioned  above  is  confirmed  by  actually  summing  the 
scattered  pulses  generated  in  the  layered  model.  Fig.  13  is  the  result.  100  scatterers  are  distri¬ 
buted  randomly  with  the  scattering  distance  between  300  km  and  400  km  in  the  CUS  model. 
A  Q  model  10U0Q  is  also  used.  The  slope  from  this  synthetic  coda  is  very  close  to  that 
predicted  the  first  trace  in  Fig.  12.  At  these  scattering  distances,  the  ratio  of  the  coda  lapse 
time  to  the  Lg  propagation  time  is  greater  than  2 


8.  Numerical  Experiment  U  Body  H  are  Scattering 

The  theory  of  body-wave  scattering  m  the  earth  has  been  widely  discussed  from  the 
points  of  view  of  scattering  from  a  discrete  volume  (Knopoff,  1959;  Wu,  1984)  or  as  superposi¬ 
tion  of  scattered  body-waves  (Sato.  1981)  These  theories,  however,  use  a  plane  wave  reacting 
with  the  scatterer  in  a  whole  space  or  .it  most  a  halfspace  In  this  section,  we  will  discuss 
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body-wave  scattering  in  the  layered  medium.  The  layered  earth  structure  is  an  important  fac¬ 
tor  which  provides  a  circumstance  for  waves  to  react  with  inhomogeneities  by  trapping  waves 
in  the  layers.  We  will  simply  extend  the  theory  developed  for  surface- waves,  i.e.,  Lg  waves,  to 
body-waves  by  using  ’locked  mode’  approximation  (Harvey,  1981).  This  is  done  by  adding  a 
rigid  cap  layer  at  great  depth  in  the  original  model  to  trap  the  body-waves.  In  doing  this,  the 
body-waves  become  a  subset  of  the  surface-wave  field.  Hence,  we  can  easily  extend  the  first- 
«  order  scattering  theory  to  cover  this  part  of  scattering  signals. 

If  care  is  taken  in  the  choice  of  the  cap  layer  properties,  the  body-wave  scattering 
waveform  will  be  generated  as  shown  in  Fig.  14.  The  second  trace  in  this  figure  is  an  exag¬ 
geration  of  the  body-wave  part  from  the  first  trace  which  is  the  scattering  wave  generated  by 
the  locked  mode  approximation  using  W'US  model.  The  third  trace  is  the  scattering  wave 
from  surface-wave  mode  only.  Body- wave  scattering  is  strong  around  the  P  and  S  arrivals. 
We  can  also  see  some  body-wave  and  surface-wave  scattering  interaction. 

From  their  spectra  plots,  we  can  see  that  the  scattered  pulses  contain  similar  frequencies 
if  there  is  no  attenuation.  However,  the  attenuation  makes  the  difference.  Fig.  14(b)  and 
14(c)  shows  the  effects  of  intrinsic  Q  and  the  scattering  Q,  respectively.  It  is  interesting  to  see 
that  the  body-wave  scattering  suffers  much  less  attenuation  than  surface-wave  scattering  espe¬ 
cially  in  high  frequency  range.  This  can  be  seen  more  apparently  in  Fig.  15  where  spectra 
under  different  Q  conditions  are  plotted  for  surface  and  body-wave  scattering  respectively  for 
a  scattering  distance  of  250  km..  The  high  frequency  spectrum  slopes  do  not  show  as  much 
variation  for  body-wave  scattering  as  for  surface-wave  scattering.  This  implies  that  the 
body-wave  scattering  retains  more  high  frequencies  than  the  surface-wave  scattering  after 
attenuation.  This  makes  the  body-wave  scattering  dominant  at  high  frequencies  as  observed 
by  Shin  (1985). 

9  Conclusion 

After  examiug  the  behavior  of  scattered  pulses  generated  in  a  layered  medium,  we  have 
the  following  conclusions: 

1)  The  ’diagonal  selection’  rule  which  predicts  no  mode  conversion  is  not  proper  for  scatter¬ 
ing  in  the  layered  elastic  medium.  Mode  conversion  must  be  considered  especially  for  the  hor¬ 
izontal  components  or  at  short  distances. 

2)  The  wave  type  .-onversion  between  Rayleigh  and  Love  waves,  on  the  other  hand,  can  be 
ignored  except  at  short  distances. 

3)  The  scatterers  near  the  surface  can  generate  well  dispersed  scattered  waves,  but  they  also 
lose  energy  faster  because  of  dissipative  attenuation  The  combined  effect  of  attenuation  and 

!  layering  makes  the  scatterer  depth  dependence  less  significant. 

4)  The  variation  of  scattering  Q  does  not  affect  the  value  of  apparent  Q  as  much  as  the 
intrinsic  Q.  For  the  same  scattering  environment,  the  attenuation  of  high  frequencies  from 
scattering  Q  does  not  seem  to  depend  much  on  the  scattering  location. 

5)  It  is  difficult  to  make  a  distinction  between  the  effect  of  scattering  Q  and  intrinsic  Q  on 
the  final  Q  value. 

6)  Larger  geometrical  spreading  decay  is  found  in  the  early  part  of  coda.  It  falls  off  by  a 
factor  of  t“°s  to  t-1  faster  than  the  later  part.  The  geometrical  correction  in  the  early 
unstable  region  of  coda  is  suggested  to  be  larger  than  the  simple  Sato’s  coircction.  This  is 
caused  mainly  by  rapid  mode  or  wave  type  conversion  from  scatterers  at  near  distances. 

7)  If  the  scattering  pulses  which  constitute  the  coda  are  dense  enough,  which  is  required  to 
make  a  smooth  decaying  coda  envelope,  the  geometrical  correction  factor  will  control  the 
value  of  coda  Q. 

8)  The  body-wave  scattering  suffers  less  attenuativc  absorption  than  the  slower  surface 
wave  scattering 

9)  The  effect  of  layering  on  the  coda  formation  is  very  important  Because  of  wave  trap¬ 
ping  in  the  layers,  even  a  weak  scatterer  distribution  may  create  observerable  coda.  This 
layering  is  the  same  reason  that  makes  the  Lg  phase  the  most  dominant  phase  in  regional 


earthquakes  seismograms. 


Acknowledgments 

This  work  was  supported  in  part  by  the  National  Science  Foundation  under  Grant  No.  CEE- 

8406577  and  by  the  AFSC  under  Contract  F19628-85-K-0029. 

References 

Aki,  K.  (1969).  Analysis  0}  the  seismic  coda  of  local  earthquakes  as  scattered  waves.  J.  Geo¬ 
phys.  Res.  7 4,  615-631. 

Aki,  K.  (1982).  Scattering  and  attenuation.  Bull.  Seism.  Soc.  Am.  72,  S319-S330, 

Aki,  K.  and  Chouet,  B.  (1975).  Origin  of  coda  waves:  source  attenuation  and  scattering 
effects.  J.  Geophys.  Res.  80,  3322-3343. 

Aki,  K.  and  Richards,  P.,  Quantitative  Seismology:  Theory  and  Method,  Vol.  1.  (\V.  H.  Free¬ 
man  and  Company,  San  Francisco  1980). 

Chernov,  L.  A.,  Wave  Propagation  in  a  Random  Medium,  (McGraw-Hill,  New  York  1960). 

Dainty,  A.  M.  (1981).  A  scattering  model  to  explain  seismic  Q  observations  in  the  lithosphere 
between  1  and  SO  Hz.  Geophys.  Res.  Letters  11,  1126-1128. 

Dainty,  A.  M.  (1984).  Influence  of  Scattering  on  Q  in  the  Lithosphere.  Final  report  to  the 
AFSOR,  Georgia  Institute  of  Technology,  Atlanta,  Georgia. 

Der,  Z.  A.,  Marshall,  M.  E.,  O’Donnell,  A.  and  McElfresh,  T.  W.  (1984).  Spatial  coherence 
structure  and  attenuation  of  the  Lg  phase,  site  effects,  and  the  interpretation  oj  the  Lg 
coda.  Bull.  Seism.  Soc.  Am.  7 4,  1125-1147. 

Harvey,  D.  J.  (1981).  Seismogram  synthesis  using  normal  mode  superposition:  locked  mode 
approximation.  Geophys.  J.  R.  Astr.  Soc.  66,  37-69. 

Haskell,  N.  A.  (1963).  Radiation  pattern  of  Rayleigh  waves  from  a  fault  of  arbitrary  dip  and 
direction  of  motion  in  a  homogeneous  medium  Bull.  Seism.  Soc.  Am.  53,  619-642. 

Haskell,  N.  A.  (1964).  Radiation  pattern  of  surface  waves  from  point  sources  in  a  multi-layered 
medium.  Bull.  Seism.  Soc.  Am.  54,  377-393. 

Hudson,  J.  (1977).  Scattering  waves  in  the  coda  of  P.  Geophys.  J.  R.  Astr.  Soc.  4 8 ,  359-374. 

Ishimaru,  A.,  Wave  Propagation  and  Scattering  in  Random  Media,  (Academic  Press,  New 
York  1977). 

Kennett,  B.  L.  N.  (1972).  Seismic  waves  in  laterally  inhomogeneous  media.  Geophys.  J.  R. 
Astr.  Soc.  27,  301-325. 

KnopofT,  L.  (1959).  Scattering  of  compressional  waves  by  spherical  obstacles.  Geophysics  24, 
30-39. 

Kopnichev,  Y.  F.  (1977).  Models  for  the  formation  of  the  coda  of  the  longitudinal  wave.  Proc. 
(Dokl.)  Acad.  Sci.  USSR  294,  13-15. 


n^wrowowo— wiwrow  wiwiwuiu  tv  vu^u«» 


-  15  - 


Malin,  P.  E.  (1978).  A  first  order  scattering  solution  for  modeling  lunar  and  terrestrial  seismic 
codas.  Ph.D.  Dissertation,  Princeton  University. 

Malin,  P.  E.  (1980).  A  first-order  scattering  solution  for  modelling  elastic  wave  codas  -  I  The 
acoustic  case.  Geophys.  J.  R.  Astr.  Soc.  6S,  361-380. 

Pulli,  J.  J.  (1984).  Attenuation  of  coda  waves  in  New  England.  Bull.  Seism  Soc.  Am.  74. 

1149-1166. 

Sato,  H.  (1977).  Energy  propagation  including  scattering  effect;  single  isotropic  scattering  J 
Phys.  Earth  25,  27-41. 

Sato,  H.  (1984).  Attenuation  and  envelope  formation  of  three-component  seismograms  of  small 
local  earthquakes  in  randomly  inhomogeneous  lithosphere.  J.  Geophys.  Res.  89.  1221- 
1241. 

Shin,  T.-C.  (1985).  Lg  and  coda  wave  studies  of  Eastern  Canada.  Ph.  D.  Dissertation,  Saint 
Louis  University. 

Wang,  C.  Y.  (1981).  Wave  theory  for  seismogram  synthesis  Ph  D.  Dissertation,  Saint  Louis 
University,  St.  Louis,  Missouri. 

Wu,  R.  S.  (1984).  Seismic  wave  scattering  and  the  small  scale  inhomogeneities  in  the  litho¬ 
sphere.  Ph.D.  Dissertation,  Massachusetts  Institute  of  Technology. 

Wu,  R.S.  and  Aki,  K.  (1985).  Scattering  characteristics  of  elastic  waves  by  an  elasteic  hetero¬ 
geneity.  Geophysics  SO,  582-595. 


Appendix 

I.  Source  Function: 

For  a  double  couple  source,  the  source  mechanism  can  be  determined  by  the  dip  d  ,  slip 
s,  and  strike  d  in  the  following  forms: 

«S  =  "k  UrR„  +  A+i  Rk  Ur)  Rdd  +  i(  Rk  U,A  Rd, 
dz  2  ci  z 

,  ,  ,  dlT,  , 

*s=  He  U#R„  -  i— -Rds 

dz 

where 

Rm  =  — sind  coss  sin2d  —  —  sin2d  sins  cos2d 

2 

Rdd  =  sins  sin2d 

Rdi  =  cos2d  sins  sind  —  cosd  coss  cosd 
R„  =  sind  coss  cos2d  —  — sin2d  sins  sin2d 
Ra.  =  cosd  coss  sind  +  cos2d  sins  cosd 


For  an  explosive  source,  only  the  Rayleigh  wave  has  a  source  function: 

dU. 


*5  =  0. 


II.  Scattering  angular  integration: 

The  angular  integrations  for  determining  the  scattering  coefficient  (Equation  7)  are 
©4  =  C0C| 

©2  =  A^af  +  C0cf  +  2C0c2c0  +  E0(a1c,+a0c2)  +  F0b0c2 
©o  =  J^oa0-  +  B0b0-  +  C0Cq  +  Duaobo  +  EoaoCo  +  F0b0c0 

where 

Ao  =  (vpw2)2  B0  =  ((2+v)X)2  CQ  =  ((2+v)/i)2 

D0  =  — 2v(2+v)pX_r  E0  =  -2v(2+v)/>/iur  F0  =  2(2+v)2X/r 

And  a’s,  b’s,  and  c’s  for  RR  conversion  are 

al  =  ^  rm^  rn  *0  FImUJ0 

bo  =  (kmUrm-c)JUI(n)  (kDt;„-0,Uim) 

Co  —  -kmErm  k0Lrn  C|  =  (kmUim+<91Urm)(knL  in+5,Uro) 
c0  =  2d,U„0IUlo 


for  RL, 


al  Dnnli/,,  ao  —  0 

bo=  o 

c2  =  lcmBmik11U<|1  C,  =  ^jUnn^jU^+kjnUjmdjU#,  C0  =  0 


for  LR, 


ai  =  D<mUfl 


ao  —  © 
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Table  1 

Central  United  States  Model  (CUS) 


d 

a 

0 

P 

Q* 

Q* 

1 

5.00 

2  89 

2.5 

600 

300 

9 

6.10 

3.52 

2.7 

600 

300 

10 

6.-10 

3.70 

2.9 

600 

300 

20 

6.70 

3.87 

3.0 

4000 

2000 

8.15 

4.70 

3.4 

4000 

2000 

Table  2 


West  United  States  Model  (WUS) 

d 

a  0  P  Qa  Q* 

2 

3.55 

2.06 

2.20 

170 

85 

3 

6.15 

3.27 

2.79 

300 

150 

18 

6.15 

3.57 

2.79 

300 

150 

8 

6.70 

3.93 

2.97 

1000 

500 

8 

6.70 

3.73 

2.97 

1000 

500 

7.80 

4  41 

3.35 

2000 

1000 

Vs< 

Ml 


Figure  1.  The  scattering  of  waves  in  the  layered  medium.  The  source,  S,  scattercr  and  receiver  R 
are  shown.  For  a  given  source-rcceivrr  distance,  the  locus  of  all  scatterers  at  a  given  depth  forms 
an  ellipse,  which  is  displayed  at  the  surface. 
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Figure  2.  A  comparison  of  waveforms  of  direct  surface  waves  and  the  scattered  waves.  The 
numbers  at  the  end  of  each  seismogram  indicates  the  source-receiver  distance  for  surface-waves  and 
the  source-  scatterer-receiver  distance  for  the  scattered  waves.  The  scattered  waveform  is  much 
more  concentrated  than  the  direct  wave 


Figure  3.  Wave  type  conversion  test.  The  upper  trace  in  each  pair  of  seismograms  does  not  include 
wave  type  conversion  and  the  lower  one  does.  The  tests  are  for  scatterers  at  the  distances  of  110 
km  and  300  km,  and  at  the  angles  0,  90  and  150  degrees  on  an  ellipse.  Note  the  waveform 
difference  in  wide  angle  scattering  (90  or  150)  for  the  short  distance  (110  km). 
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Figure  4.  Mode  conversion  test.  Cases  discussed  are  for  mode  conversion  permitted  for  up  to  1,  5 
and  10  neighboring  modes.  Three  scattering  distances  110  km.  300  km  and  500  km  are  included. 
Two  columns  correspond  to  cases  with  or  without  wave  type  conversion  respectively  Plots  a,  b 
and  c  are  for  three  different  components  for  the  horizontal  components.  Obvious  mode  conversion 
occurs  for  the  short  distance  scattering  (110  km). 
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Figure  4c. 
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Figure  6.  A  comparison  of  the  effect  of  scattering  Q  and  intrinsic  Q.  The  scattered  are  at  the 
depth  of  0.5  km  and  the  scattering  distance  is  300  km.  101 ‘JQ  means  10  scattering  attenuation 
unit  (l!,»)  and  double  anelastic  attenuation  as  given  by  Q  model  in  Table  1  Note  that  the  fre¬ 
quency  content  of  the  nUOQ  cases  is  very  similar.  The  effect  of  scattering  Q  seem-  to  enhance  the 
high  frequency  and  higher  mode  signal- 
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Figure  7.  The  scattering  waveforms  and  their  spectra  under  different  attenuation  conditions.  The 
scattering  distance  is  500  km  and  the  scatterer  is  at  the  depth  of  5  km  (a)  and  15  km  (b)  in  Cl’S 
model.  The  spectra  plotted  are  in  the  same  sequence  as  the  seismograms  but  shifted  by  a  factor  of 
ten  for  clarity.  Their  high  frequency  spectra  slopes  are  also  indicated.  After  traveling  large  dis¬ 
tances,  the  high  frequencies  are  absorbed  due  to  scattering  and  anelastic  absorption.  The  scatter¬ 
ing  Q  seems  to  have  similar  effect  for  scatterers  from  different  depths,  however,  intrinsic  Q  affects 
the  shallower  scatterer  more  since  the  intrinsic  Q  value  is  low  near  the  surface.  In  case  (a),  the 
anelastic  Q  controls  the  decay  of  wave  energy.  On  the  other  hand,  in  case  (b)  the  scattering  Q  is 
dominant  in  controlling  the  frequency  decay. 
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Figures  8.  The  geometrical  spreading  factor  for  scattering  waves  in  the  layered  medium  Displays 
are  for  Z  component  (a)  and  SH  component  (b)  with  epicentral  distance  at  100  km  in  Cl'S  model 
The  scattering  distance  varies  from  110  km  to  GOO  km  The  top  curve  is  from  the  average  peak- 
to-peak  amplitude  of  composed  scattering  waveforms,  and  the  lower  two  curves  are  for  amplitudes 
after  narrow  bandpass  filtering  at  at  1  Mr  and  3  Hz  respectively  No  Q  has  been  applied  Because 
of  mode  conversion,  we  can  see  much  faster  amplitude  decay  in  the  front  part  (<  200km)  The 
slopes  of  piecewise  linear  segments  fit  to  the  data  are  shown 
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Figure  9.  Scattering  geometrical  spreading  in  the  Wl'S  model.  Two  components,  vertical  and 
tangential,  are  displayed,  (a)  is  for  source  -receiver  distance  at  300  km  and  (b)  at  100  km  Note 
that  the  slope  changes  at  600  km  for  (a)  and  at  200  km  for  (b)  The  WL'S  model  has  more  rapid 
geometrical  decay  rate  than  the  Cl'S  model.  The  slope  difference  between  the  two  regions  is  also 

larger. 
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Figure  10  A  modification  of  geometrical  spreading  correction  of  Sato’s  surface-wave  scattering 
model  A  higher  decay  rate  of  t-04  to  t-1  due  to  mode  conversion  and  wave-type  conversion  is 
assumed  for  the  lapse  time  ratio  smaller  than  2. 
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Figure  11.  Duration  of  scattering  pulse  at  different  distances.  A  Q  model  10U1Q  in  CIS  model 
has  been  applied.  The  upper  and  lower  sets  of  data  are  from  SH  and  Z  synthetics,  respectively 
An  average  duration  expansion  rate  0.05  with  distance  is  obtained. 
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'igure  12.  Superposition  of  scattering  pulses  to  form  the  coda  envelope.  The  part  (a)  simulates  the 
astern  United  States  case  Spreading  factor  m  largely  affects  the  coda  decay.  Also  note  that  the 
nalysis  of  first  part  (I)  and  the  rest  (D)  yields  different  apparent  Q  values  The  part  (b)  simulates 
he  western  United  States  case.  Since  the  pulse  superposition  has  been  saturated  by  the  rapid 
ecay  due  to  Q,  the  inferred  Q  is  not  dependent  upon  the  shape  factors  of  the  individual  scattered 
ulses  much. 
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Figure  13.  Coda  synthetics  after  summing  over  100  scattered  pulses  which  are  from  scatterers  dis¬ 
tributed  evenly  between  scatterer  distances  300  km  and  400  km.  A  Q  model  10U1Q  has  been 
applied.  The  decay  of  the  coda  follows  the  shape  predicted  by  the  pulse  summation. 
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Figure  14.  The  spectra  of  body-wave  and  surface-wave  scattered  pulses,  (a),  (b),  and  (c)  are  for 
three  different  attenuation  conditions.  Three  traces  in  each  plot  are  for  complete,  body-wave  only, 
and  surface-wave  only.  The  spectra  of  surface-wave  only  is  shifted  down  by  a  factor  of  100  for 
clarity.  Note  that  the  frequency  content  is  similar  (see  (a))  for  body-  or  surface-wave  scattering  if 
no  Q  effect  is  included,  (b)  shows  the  case  under  scattering  attenuation.  The  attenuation  of 
body-wave  scattered  waveforms  is  less  than  the  surface-wave  scattered  arrivals,  (c)  is  for  the 
intrinsic  Q.  Again,  the  body  wave  scattering  shows  stronger  resistance  to  attenuation. 
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Figure  15.  The  decay  of  high  frequency  component  under  different  attenuation  conditions  for 
surface-wave  and  body-wave  scattering  fcr  a  scattering  distance  of  250  km..  The  body-wave 
scattering  shows  much  smaller  attenuation  compared  to  surface-wave  scattering  attenuation.  This 
may  explain  much  higher  frequencies  are  observed  from  P-coda  or  S-coda  Each  spectrum  has  been 
shifted  by  a  factor  of  ten  for  clarity. 
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APPLICATION  OF  FREQUENCY  VARIABLE  FILTERS 
TO  SURFACE- WAVE  AMPLITUDE  ANALYSIS 

By  David  R.  Russell,  Robert  B.  Herrmann,  and  Horng-Jye  Hwang 


ABSTRACT 

The  problem  of  spectral  biasing  due  to  frequency  domain  filtering  of  surface-wave 
seismograms  is  investigated,  and  the  method  of  frequency  variable  filters  (FVF)  is  developed 
t>  compensate  for  this  bias.  As  a  result,  the  FVF  can  significantly  improve  signal  to  noise  in 
t  ie  filtering  process,  except  at  points  which  require  increased  frequency  domain  resolution  due 
to  biasing.  A  detailed  comparison  of  some  currently  accepted  surface-wave  filters  is  made  in 
order  to  clarify  the  development  of  the  FVF  algorithm.  Three  long  period  surface-wave 
seismograms  are  tested  with  FVF  and  compa-ed  to  two  other  methods,  the  multiple  filter 
technique  (MFT)  and  the  phase  matched  filter  (PMF).  Emphasis  is  placed  on  finding  limita¬ 
tions  in  all  the  methods,  not  on  routine  processing.  Results  of  the  tests  show  that  the  FVF 
and  PMF  are  an  improvement  over  MFT,  in  that  the  results  of  processing  can  be  diagnosed  in 
both  the  time  and  frequency  domains.  In  addition,  the  FVF  is  more  successful  than  PMF  in 
removing  higher  mode  contamination  from  fundamental  mode  Rayleigh  waves.  Another  con- 
clusion  is  that  care  must  be  taken  in  defining  the  limits  of  frequency  domain  resolution  for  the 
FVF.  Some  bias  due  to  convolutional  smoothing  must  be  accepted,  or  the  FVF  becomes  a 
simple  all-pass  filter.  MFT  and  PMF  have  a  similar  problem,  which  is  due  fundamentally  to 
the  uncertainty  principle  between  time  and  frequency. 

Introduction 

Recently ,  attention  has  been  focused  on  the  problem  of  filtering  seismograms  to  isolate 
propagating  normal  modes,  especially  in  the  context  of  determining  spectral  amplitudes  for 
nuclear  yield  estimation.  Using  the  multiple  filter  technique  (MFT)  (Dziewonski  et  al. 
1969),  Yacvub  (1983)  determined  average  spectral  amplitudes  around  the  20  second  period 
Herrin  and  Goforth  (1977)  developed  the  phase-matched  filter  (PMF)  to  refine  group  and 
phase  velocities  of  normal  modes,  and  Stevens  (1986a)  used  this  method  to  isolate  fundamen¬ 
tal  mode  spectra  across  the  entire  observed  bandwidth.  Hwang  and  Mitchell  (1986)  used  the 
time-variable  filter  (TVF)  developed  by  Landisman  et  al.  (1969)  to  isolate  spectral  ampli¬ 
tudes  for  inter-station  Green’s  functions.  Russell  et  al.  (1986)  combined  the  PMF  and  TVF 
methods  to  form  a  frequency-variable  filter  (FVF),  in  order  to  improve  spectral  amplitude 
estimates  of  explosion  generated  waveforms. 

This  paper  will  extend  the  FVT  method  in  order  to  reduce  problems  of  spectral  biasing, 
and  will  discuss  in  some  detail  relationships  between  the  above  filters. 

General  theory 

For  surface-wave  studies,  the  objective  is  to  process  s(t).  the  raw  surface-wave  time  his¬ 
tory,  by  a  series  of  convolutional  filters  designed  for  modal  isolation.  In  the  context  of  this 
discussion,  the  definition  of  "convolution"  is  extended  to  include  the  general  case  of  frequency 
varying  filters,  which  have  a  changing  bandwidth  during  the  convolution  operation.  This  will 
produce  a  time  trace  (or  frequency  spectrum)  which  can  be  used  for  further  analysis. 

Most  surface-wave  filters  can  be  expressed  by  the  following  relations: 

00 

*  to  =  f  PjM  W^t.w)  S(w)  ef^dw  (1) 
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MO  =  ^  /  Pf ‘M  e'w'd^  (3) 

*  -oc 

where  ty(ut)  =  Fourier  transform  of  V’  (t)  , 

SM  =  £  IS^w)  |  e'‘kjM‘  +  N(-') 

j 

P,(w)  =  eU‘)H‘  ,  P-‘M  =  e^1* 


S(u>)  is  the  total  spectrum  of  the  seismogram,  composed  of  a  sum  of  normal  modes  and 
an  additive  component  N(-j),  corresponding  to  multi-pathed  signals,  interfering  events,  and 
incoherent  noise.  The  purpose  of  the  above  filters  is  to  isolate  the  j’th  mode  of  interest  so  the 
amplitude  spectrum  |Sj(u)  |  and  wavenumber  spectrum  k,(w)  can  be  recovered,  with  the 
understanding  that  the  spectra  may  contain  both  source  and  path  contributions.  This  is 
probably  too  idealistic  a  goal  in  some  cases,  as  pointed  out  by  Der  (1986).  Due  to  scattering 
and  reflections,  there  may  not  be  a  pure  isolated  mode  to  recover.  A  more  accurate  statement 
is  that  the  purpose  of  such  filters  is  to  isolate  seismic  energy  propagating  in  the  vicinity  of 
desired  modes  of  interest. 

Pj(w)  is  a  phase-matched  filter.  The  wavenumber  estimate  kj  should  be  near  the  true 
modal  wavenumber,  in  order  to  compress  the  energy  of  the  desired  signal  about  zero-lag  in  the 
time  domain,  forming  a  "pseudo-autocorrelation  function"  rp  (t).  Herrin  and  Goforth  (1977) 
discussed  this  in  detail. 

Wf(t,w)  is  a  time  and  frequency  variable  window  used  to  isolate  modes  of  interest  and 
improve  signal  to  noise  ratios.  It  is  symmetric  about  position  r  in  the  time  domain,  with  a 
width  controlled  by  the  frequency  w. 

H(u>— w0)  is  a  frequency  domain  convolution  filter,  used  to  isolate  the  energy  of  pro¬ 
pagating  normal  modes  in  multiple  filter  analysis. 


Review  of  surface-wave  applications 

Various  combinations  of  the  above  filters  have  been  used  for  modal  isolation,  four  of 
which  will  be  discussed  in  detail  below. 


Case  1:  P(w)  m  1  ,  \V,(t,w)  =  1  .  This  is  the  basis  for  the  multiple  filter  technique 

(MFT)  (Dziewonski  et  al.  ,  1969).  In  this  case  ^(u;)  =  S(w),  so  the  raw  spectrum  S(— •)  is 
input  into  (2),  which  is  the  MFT  evaluated  at  u»0.  H(w— w0)  is  a  narrow  bandpass  filter  (usu¬ 
ally  a  band-limited  Gaussian),  symmetric  about  ~0.  The  symmetry  about  positive  frequencies 
causes  the  time  signal  h(t)  to  be  complex,  with  the  modulus  having  maxima  at  the  group 
velocities  of  the  signal  modes.  Herrmann  (1973)  showed  that  under  the  condition  of  an 
approximately  flat  amplitude  spectrum  and  linear  phase  delay  of  the  j’th  mode  across  the 
width  of  H(w-u0), 


2>ra  >/2 


E  |Si(w'o)  I  exp[i(w0t-kj(u;0)x)]  exp  - 


where  o  is  a  constant  controlling  the  width  of  the  Gaussian  filter,  and  t,  is  the  group  delay  of 
the  j’th  mode.  Evaluating  (4)  at  multiple  frequencies  will  extract  the  spectrum  of  the  j’th 
mode,  if  it  is  suitably  smooth. 


Case  2:  \V,(t,uj)  =  1  .  This  is  the  basis  for  the  "residual  dispersion"  technique  of 

Dziewonski  et  al.  (1972).  The  MFT  is  first  applied  to  determine  wavenumber  estimates 
from  the  instantaneous  phases  of  (4),  in  order  to  construct  the  phase-matched  filter  P(-0 
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Equation  (1)  now  represents  the  seismogram  with  the  dispersion  removed  from  the  j'th  mode. 
The  signal  is  Fourier  transformed  and  the  MFT  applied  again  (2).  This  method  was  intro¬ 
duced  to  remove  biasing  caused  by  the  signal  phase  in  MFT  due  to  dispersion. 

Case  3:  P(w)  =  1,  Wr(t,w)  =  W^jftjU/)  .  This  is  the  time-variable  filter  due  to 

Landisman  et  al.  (1969).  MFT  is  performed  first  to  determine  group  delays  tj(u.’)  ,  and  the 
window  \VT  is  symmetrically  centered  about  the  group  delay  in  equation  (1).  The  width  of 
the  window  is  a  multiple  of  the  period  of  interest.  V'  (t)  should  therefore  contain  only  energv 
associated  with  the  mode  of  interest,  and  so  (1)  will  be  equivalent  to  the  time  domain  signal 
of  the  isolated  mode. 


Case  4:  W,(t,u)  —  Wq(i)  .  This  is  the  phase-matched  filter  metliod  of  Herrin  and 

Goforth  (1977).  MFT  is  performed  first  to  estimate  wavenumbers  for  the  phase-matched  filter 
P(w).  This  differs  from  the  residual  dispersion  technique  (case  2),  in  that  the  phase  of  the  sig¬ 
nal  is  found  by  integrating  the  group  delay  calculated  in  MFT.  The  window  W^t)  is  no 
longer  frequency  dependent,  so  it  can  be  factored  from  the  integral  in  (1).  In  the  time  domain 
W„(t)  is  centered  about  zero-lag,  and  should  isolate  the  energy  of  the  desired  mode  of  interest. 
In  the  frequency  domain,  (1)  will  be  a  true  convolution  of  the  pseudo-autocorrelation  spec¬ 
trum  with  the  frequency  equivalent  of  W0(t).  Equation  (3)  is  the  time  domain  signal  of  the 
isolated  mode. 

All  of  the  above  methods  attempt  to  isolate  spectral  modes  of  interest  via  frequency 
domain  convolutions.  As  Dziewonski  and  Hales  (1972)  pointed  out,  the  process  of  convolution 
will  distort  the  signal  spectrum  unless  the  prescribed  filters  are  Dirac  impulses  in  the  fre¬ 
quency  domain.  The  process  of  residual  dispersion  (case  2),  or  phase-matched  filtering  (case  4) 
can  minimize  phase  distortion,  but  they  do  not  address  the  problem  of  amplitude  distortion, 
or  bias  due  to  the  convolution  operation. 

FVF  theory 

Let  W^t,w)  =  Wo(t,w).  This  defines  the  method  of  frequency  variable  filtering  (FVF).  It 
combines  the  beneficial  aspects  of  a  time-variable  filter  (TVF)  and  a  phase-matched  filter 
(PMF),  and  allows  the  analyst  to  place  maximum  error  bounds  on  the  amount  of  tolerable 
amplitude  bias.  The  raw  spectrum  S(ui)  is  phase-matched  filtered  to  compress  the  energy  of 
the  signal  of  interest  about  zero-lag  in  the  time  domain.  Then,  each  harmonic  component  of 
this  pseudo-autocorrelation  function  is  windowed  about  zero-lag  with  W0(t,w),  with  the  width 
of  the  window  proportional  to  the  period  of  the  individual  harmonic.  These  operations  are 
done  in  equation  (1),  which  is  then  Fourier  transformed  for  the  phase-matched  spectrum  of 
interest.  Equation  (3)  replaces  the  phase  removed  in  (I)  to  give  the  time  domain  signal  of  the 
isolated  mode.  The  technique  is  quite  similar  to  the  time-variable  filter  (case  3),  except  that 
the  phase  of  the  mode  of  interest  is  first  removed  to  minimize  spectral  distortions  due  to 
phase  fluctuations.  It  should  be  noted  that  in  the  context  of  residual  dispersion  measure¬ 
ments,  this  method  was  first  recommended  by  Dziewonski  and  Hales  (1972). 

An  advantage  of  FVF  is  that  the  time  window  \V0(t,w)  can  be  constructed  at  each  fre¬ 
quency  to  control  the  amount  of  bias  due  to  windowing.  The  exact  form  of  the  window  will 
be  defined  below,  but  first  a  discussion  of  the  problem  of  biasing  is  presented. 

In  calculating  equation  (1),  it  is  assumed  that  the  window  W^t.w)  does  not  distort  the 
spectrum  of  the  j’th  mode.  This  is  not  true  in  practice,  and  the  purpose  of  this  section  is  to 
approximately  calculate  the  bias  in  the  frequency  domain  due  to  time  domain  windowing. 
The  analysis  is  similar  to  Jenkins  and  Watts  (1968,  p.  247).  From  equation  (1),  define  the 
signal  pseudo-autocorrelation  function  as  a  phase-matched  single-mode  signal  uncontaminated 
by  noise: 

00 

V-j(t)  =  ^  /  ISM  |  ei&M  *  e“iw'dw  (5) 
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where 

*M  =  kj(^)  -  kj(w) 

is  the  residual  wavenumber  left  from  the  phase-matched  filter  process.  The  bias  in  the  fre¬ 
quency  domain  will  be  the  transformed  difference  between  the  windowed  signal  pseudo¬ 
autocorrelation  function  and  (5): 

OO 

BM  =  /  (\v0(t)  -  1]  ^(t)  e-^'dt.  (6) 

— OO 


Two  windows  are  examined  for  bias:  the  cosine  and  Parzen  (also  called  de  la  Valle  - 
Poisson)  windows.  The  cosine  window  is  defined  as: 


w0c(t) 


cos(rt/2T)  |t|<T 
0  |t  |  >T 


(7) 


where  T  is  the  one-sided  width  of  the  cosine  window.  The  maximum  sidelobe  level  for  the 
cosine  window  is  -23  dB  (7  percent)  of  the  main  lobe  (Harris,  1978).  The  Parzen  window  is 
defined  as: 


W0p(t)  = 


-  6(t/T)2  +  B(  |t  |/T)3  |t|<T/2 


2(1  -  1 1 1 / T)3  T/2<  |t|<T 


0 


It  |  >T 


(8) 


where  T  is  the  one-sided  width  of  the  Parzen  window.  The  maximum  sidelobe  level  for  the 
Parzen  window  is  -53  dB  (0.2  percent)  of  the  main  lobe  (Harris,  1978).  For  a  given  width  T, 
the  Parzen  window  is  a  much  smoother  convolutional  filter  in  the  frequency  domain  due  to 
low  sidelobes. 

To  calculate  the  approximate  bias,  substitute  the  windows  into  (6)  and  keep  terms  only 
on  the  order  of  1/T2  or  more.  For  the  cosine  window, 


OO 


Bc(w) 


ir 

8T2 


/  -t2  0,(t)  e-**  dt  +  0(1/T4) 


and  the  Parzen  window, 

OO 

BPM  =  /  'l2  W  e--‘  dt  +  0(1/T3) 

^  — OO 


(9) 


Notice  that  the  infinite  limits  are  kept  for  the  integrals.  This  is  under  the  assumption  that  T 
is  wide  enough  to  insure  an  insignificant  truncation  of  the  signal  pseudo-autocorrelation  func¬ 
tion,  resulting  in  a  negligible  signal  outside  the  limits  +T,  -T  .  Making  use  of  Fourier 
transform  properties  of  differentiation  (Papoulis,  1962,  p.  16)  gives 

B'M  =  —  */'(-) ,  B’M  =  */'M  (10) 

where  the  double  primes  indicate  the  second  derivative  with  respect  to  angular  frequency  of 
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the  signal  pseudo-autocorrelation  spectrum  Taking  the  ratio  of  Bp  to  Bc  shows  that  the  Par- 
zen  window  has  almost  five  times  the  bias  of  the  cosine  window,  which  suggests  that  it  may 
be  a  poor  choice  for  a  windowing  function.  However,  further  investigation  of  the  spectral 
second  derivative  yields  interesting  results. 

From  equation  (5),  the  spectrum  of  the  signal  pseudo-autocorrelation  function  is 

W  =  |Sj(w)  I  ef*M«  =  oei#.  (11) 


The  terms  a  and  0  are  chosen  for  notational  convenience, 
of  the  signal  pseudo-autocorrelation  spectrum  gives 


♦/'(<*»)=  [a”  -a(/)V  +  [at"  +2aV]ei<,  +  ,"2>  . 


Calculating  the  second  derivative 


(12) 


If  the  phase-matched  filter  is  successful,  the  difference  between  the  estimate  and  true 
wavenumber  ^  )  should  approach  zero.  Therefore,  keeping  only  first  order  residual  phase 

terms  for  9,  6  ,  and  0  ,  the  Parzen  window  amplitude  bias  is 


and  the  phase  bias  is 


Bp<  = 


6(a6"  +  2aV) 

T2a  +  6a" 


(13) 


(14) 


A  similar  expression  for  the  cosine  window  amplitude  bias  is 


and  for  phase  bias, 


_2 

Bc“  =  -—z-a" 
8T2 


Dti  =  ^(aO"  +  2a  V) 

8T2o  +  jrV' 


(15) 


(16) 


Equation  (14)  and  (16)  show  that  the  bias  in  phase  is  independent  of  constant  phase 
values,  and  that  if  the  first  and  second  derivatives  are  small,  the  phase  bias  will  also  be  small. 
Therefore,  it  is  advisable  to  use  smooth  (low  sidelobe)  windows  such  as  the  Parzen  for  calcu¬ 
lating  residual  phases  in  the  matched  filtering  process.  However,  it  should  be  clear  that  spec¬ 
tral  amplitudes  may  be  quite  biased  when  there  is  significant  curvature  in  the  amplitude  spec¬ 
trum,  corresponding  to  a  large  spectral  second  derivative,  a  '  .  Examples  of  this  are  narrow¬ 
band  spectra  and  the  vicinity  of  sharply  changing  band  edges.  The  fact  that  phase-matched 
filtering  is  an  iterative  process  can  reduce  residual  bias  in  phase,  but  this  will  not  help  the 
bias  in  the  amplitude  spectrum,  since  it  is  to  first  order  independent  of  phase. 

To  implement  the  FVF  algorithm,  first  note  that  the  4'j"  term  in  equation  (10) 
corresponds  to  the  true  signal  pseudo-autocorrelation  function,  which  is  unknown.  However, 
with  the  same  analysis  used  to  calculate  (9),  it  can  be  directly  shown  that 

B‘M=  +0(1/T<)  (17) 


where  the  superscript  ”  p  ”  indicates  the  second  derivative  of  the  Parzen  windowed  pseudo- 
autocorrelation  function  defined  by  (1).  This  is  under  the  assumption  that  the  effect  of  noise 
on  the  windowed  function  is  negligible,  which  is  not  unreasonable  since  the  Parzen  window 
width  used  for  phase  processing  can  be  relatively  narrow  (see  (14)  above).  However,  care 
should  be  taken  in  using  the  algorithm  for  seismograms  with  extremely  high  random  noise 
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levels. 


Define  the  maximum  tolerable  bias  error  relative  to  the  maximum  signal  amplitude  as 

|B|m„ 


F  = 

‘-'mu 


I  I  mu 


For  the  cosine  window,  substitute  (18)  into  (17)  and  solve  for  T,  the  one-sided  width  of  the 
window: 


TM  = 


T2  2 

^mu  I  'J'/'  I  mu 


Notice  that  the  width  is  now  frequency  dependent.  Equation  (19)  defines  the  "curvature 
correction"  for  reducing  bias  in  the  FVF.  To  account  for  zero-crossings  of  the  second  deriva¬ 
tive,  let 

2  ttC 

Tc(w)  =  max  ,  T(u)  (20) 

tjj 


where  C  is  a  constant  multiple  of  the  period  of  interest.  Using  equations  (7)  and  (20),  the 
instantaneous  cosine  window  for  the  FVF  algorithm  is  defined  as 


W0'(t,w)  = 


cos  [n/2TcHj  |t  |<rM 

0  |t|>rM. 


The  FVF  algorithm  can  now  be  formally  developed  as  follows. 

FVF  algorithm 

Step  1.  Using  Herrin  and  Goforth’s  iterative  phase-matched  filter  process  (1977),  isolate  the 
mode  of  interest  with  a  Parzen  windowed  pseudo- autocorrelation  function 

00 

W)  =  ~  /  PjM  w;(t)  S(ui)  e'wld-J  (22) 


and  save  the  final  wavenumber  estimate  used  in  the  phase-matched  filter 


PjM  =  ***** 


Step  2.  From  equation  (22),  calculate  the  spectrum  of  the  windowed  mode 


«7(w)~  /  0*(t)  e-^'dt 


♦'"(w)  =  /  -r  V/(t)  e— ’ ‘tit 


and  its  second  derivative 
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Substitute  equations  (24)  and  (25)  into  (19)  for  the  instantaneous  window  width  Tc(;j). 

Step  3.  Calculate  the  corrected  pseudo-autocorrelation  function  and  its  spectrum,  using  the 
phase-matched  filter  v_3)  and  the  instantaneous  cosine  window  (21): 

00 

M  /  Pj(")  W'(t.)  S (u»)  e-'dw  (26) 


oo 

W  =  /  W  e~*"ldw 

— OO 


(27) 


Step  4.  Using  the  phase-matched  filter  (23)  and  the  corrected  pseudo-autocorrelation  spec¬ 
trum  (27),  calculate  the  corrected  spectrum  of  the  isolated  mode: 


S,M  -*,H  .■***' 


(28) 


Taking  the  inverse  Fourier  transform  of  (28)  gives  the  isolated  mode  in  the  time  domain,  as  in 

equation  (3). 


Examples  and  discussion 

Three  seismic  events  were  chosen  to  illustrate  various  aspects  of  the  F\T  algorithm  and 
its  relation  to  other  convolutional  type  filters.  For  comparison,  MFT  and  PMF  were  also 
applied  to  the  events.  The  MFT  (3)  was  constructed  with  a  Gaussian  filter 


H(«)  = 


exp(-a  w  2/udi)  |w  |  <  u)c 
0  |w  |  >  wc 


The  Gaussian  cutoff  frequency  <jc  and  the  filter  width  parameter  a  were  constructed  with 
values  found  in  Dziewonski  et  al.  (1969)  and  Herrmann  (1973), 

wc  =  w0/4  a  =  16  zr  . 


The  residual  dispersion  method  (Case  2)  was  not  incorporated  into  the  MFT,  in  order  to  illus¬ 
trate  the  effects  of  phase  distortion. 

The  PMF  was  implemented  using  the  algorithm  developed  by  Herrin  and  Goforth 
(1977).  The  Parzen  window  (8)  was  used  for  phase  processing,  and  the  cosine  window  (7)  was 
used  to  extract  the  amplitude  spectrum.  The  one-sided  window  width  was  set  to  a  default 
value  of  1.5  times  the  maximum  period  found  by  NfFT  analysis. 

The  FVF  was  constructed  using  the  algorithm  developed  above.  Parameters  EmlJ  and  C 
in  (18)  and  (20)  were  set  to 

Em„  =  .035  C  =  2.5  . 

This  defines  a  minimum  one-sided  cosine  window  width  of  2.5  times  the  period  of  interest, 
and  a  maximum  bias  error  of  3.5%  .  As  stated  above,  some  error  must  be  tolerated  in  the 
filtering  process.  Setting  Em„  =  0  would  simply  turn  the  F\T  into  an  allpass  filter.  The 
problem  of  bias  error  is  also  inherent  in  both  PMF  and  MFT  filters,  as  will  be  shown.  The 
advantage  of  FVF  is  that  the  bias  can  be  controlled. 
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The  sample  seismograms  were  chosen  to  illustrate  and  contrast  the  above  filters,  and 
should  not  be  considered  as  typical"  events  for  processing.  The  seismograms  were  generated 
by  earthquake  sources,  and  no  attempt  was  made  to  remove  the  spectral  effects  of  the  source 
mechanism  or  the  source  time  function  from  the  events.  However,  instrument  deconvolution 
was  performed. 

Event  I  is  a  synthetic  seismogram  recorded  at  a  distance  of  1000  kilometers  from  a  45 
degree  pure  dip-slip  source  mechanism  at  a  depth  of  50  kilometers.  The  seismogram  is  com¬ 
posed  of  vertical  component  fundamental  and  first  higher  mode  Rayleigh  waves.  Figure  1 
(top)  is  the  amplitude  spectrum  of  the  signal,  showing  clear  spectral  contamination  of  the  fun¬ 
damental  mode  by  the  first  higher  mode  at  periods  less  than  20  seconds.  The  bottom  figure  is 
the  multiple  filter  contour  map  (Dziewonski  et  al  ,  1969)  giving  the  distribution  of  seismic 
energy  as  a  function  of  period  and  group  velocity. 

For  MFT  processing,  maximum  values  on  the  contour  map  were  picked  for  the  funda¬ 
mental  mode,  and  spectral  amplitudes  were  calculated  using  equation  (4).  Group  velocities 
were  picked  from  the  maximum  amplitude  points  on  the  contour  map  and  used  as  input  for 
PMF  and  FVF  processing. 

Figure  2(a)  shows  the  final  pseudo-autocorrelation  functions  for  PMF  and  FVF  process¬ 
ing.  A  cosine  window  with  a  one-sided  width  of  22  seconds  was  required  to  exclude  the  higher 
mode  in  the  PMF,  and  this  width  was  also  used  for  processing  bias  estimates  in  FVF.  Notice 
that  the  width  of  the  FVF  pseudo-autocorrelation  function  appears  to  be  slightly  wider  than 
the  PMF  in  Figure  2(a).  Figure  12  shows  the  actual  window  width  as  a  function  of  period. 
The  FVF  window  changes  as  a  multiple  of  the  period  of  interest  (20),  and  in  this  case,  no  cur¬ 
vature  correction  (19)  was  necessary. 

Figure  2(b)  shows  the  time  domain  results  of  PMF  and  FVF,  respectively.  The  isolated 
fundamental  mode  for  each  process  appears  almost  identical.  Figure  3  is  more  informative, 
contrasting  the  difference  between  the  theoretical  fundamental  mode  and  the  filtered  ampli¬ 
tude  spectra  for  the  three  processes.  The  MFT  and  FVF  are  almost  identical  to  the  theoreti¬ 
cal  spectrum,  while  PMF  exhibits  biasing  over  the  entire  spectrum.  It  is  most  pronounced  at 
30  seconds,  which  corresponds  to  the  region  of  highest  curvature. 

Event  I  is  somewhat  extreme,  in  that  the  presence  of  the  higher  mode  forces  a  narrow 
bandwidth  for  the  PMF  window,  causing  a  large  bias  as  predicted  by  equation  (10).  In  prac¬ 
tice,  when  dealing  with  smoothly  varying  amplitude  spectra  due  to  shallow  depth  explosions, 
higher  modes  may  not  be  significantly  excited,  and  the  PMF  can  be  constructed  with  time 
windows  exceeding  100  seconds.  This  was  graphically  demonstrated  in  Stevens  (1986b).  Long 
period  spectral  curvature  was  slight  for  the  event  Stevens  chose,  and  there  was  no  apparent 
long  period  bias  for  two-sided  window  widths  varying  from  50  to  1000  seconds.  This  indi¬ 
cates  that  in  that  case,  PMF  was  justified  for  spectral  amplitude  determination. 

To  demonstrate  the  effect  of  sharply  varying  curvature  on  the  filters,  events  II  and  III, 
with  pronounced  spectral  nulls  due  to  double-couple  source  mechanisms,  were  picked.  Both 
seismograms  are  actual  events  recorded  on  long  period  WWSSN  vertical  component  seismome¬ 
ters.  Event  II  was  recorded  at  station  SHI  in  Iran,  from  a  source  located  in  the  Gulf  of  Aden, 
with  the  signal  propagating  across  the  southeastern  Arabian  Plate. 

This  seismogram  is  composed  of  two  superimposed  events,  the  initial  earthquak  and  a 
stronger  event  occurring  approximately  175  seconds  after  the  original.  The  multiple  filter 
contour  map  in  Figure  4  separates  the  superimposed  events  into  two  distinct  energy  bands, 
which  are  almost  identical  in  shape,  indicating  a  similar  source  mechanism  and  propagation 
path. 

Pseudo-autocorrelation  functions  for  PMF  and  FVF  analysis  are  given  in  Figure  5.  The 
windowed  results  are  almost  identical,  so  only  the  FVF  isolated  mode  is  shown  in  Figure  6. 
This  figure  graphically  demonstrates  the  power  of  matched  filters  to  effectively  separate  super¬ 
imposed  signals,  in  addition  to  improving  the  signal  to  noise  ratio. 
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Since  there  is  no  theoretical  reference  spectrum,  the  filtered  amplitude  spectra  were  sim¬ 
ply  superimposed  on  each  other  in  Figure  7  for  comparison.  The  most  significant  difference  is 
in  the  MFT  spectrum,  at  9  seconds  and  between  30-40  seconds.  At  9  seconds,  the  Gaussian 
filter  (29)  has  approximate  half-amplitude  values  (  e~7  )  at  8  and  10.2  seconds.  From  Figure 
7,  it  is  clear  that  the  spectrum  has  approximately  the  same  bandwidth  as  the  Gaussian  filter 
at  this  point,  which  violates  the  requirement  for  a  flat  amplitude  spectrum  across  the 
passband  in  equation  (4).  This  results  in  the  bserved  low  amplitude.  Between  30  and  40 
seconds,  applying  Dziewonski’s  residual  dispersio*.  method  (case  2)  causes  the  discrepancy  in 
this  period  range  to  disappear,  indicating  that  there  is  enough  slope  in  the  group  velocity  (see 
Figure  4)  to  introduce  significant  phase  distortion  into  MFT  amplitude  analysis. 

The  effect  of  curvature  corrections  (19)  in  the  FVF  can  be  seen  in  Figure  12.  There  is 
an  increase  in  window  width  at  9,  12,  and  17  seconds,  corresponding  to  the  spectral  peak'  and 
nulls  of  Figure  7.  It  should  be  noted  that  the  FVF  behaves  as  poorly  as  the  MFT  at  9 
seconds,  if  curvature  corrections  are  not  incorporated. 

Event  III  was  recorded  at  U.3.  west  coast  station  COR,  from  an  earthquake  occurring  in 
the  Aleutian  Islands.  The  distance  between  source  and  receiver  is  4005  kilometers,  and  the 
propagation  path  is  oceanic,  as  seen  on  the  MFT  group  velocity  contour  plot  (Figure  10). 
The  spectral  plot  in  Figure  8  exhibits  two  distinct  nulls  in  the  signal  passband,  at  18  and  25 
seconds.  The  null  at  18  seconds  results  in  a  strong,  narrowband  spectral  peak  at  17  seconds, 
shown  by  the  arrow  on  the  spectral  plot  in  Figure  8.  This  occurs  in  the  period  range  where 
the  group  velocity  slope  is  almost  vertical,  resulting  in  a  strongly  dispersed,  sinusoidal  time 
domain  signal.  This  can  be  seen  on  the  MFT  contour  map  by  comparing  the  energy  of  the  sig¬ 
nal  at  17  seconds  with  the  corresponding  time  domain  signal  on  the  right. 

The  event  was  picked  to  illustrate  how  all  of  the  above  filtering  methods  can  fail,  if  the 
entire  passband  is  considered  to  be  the  desired  signal.  This  is  demonstrated  by  the  PMF  and 
FVF  pseudo-autocorrelation  functions  in  Figure  9.  The  windowed  functions  delete  a  consider¬ 
able  portion  of  the  signal  energy,  seen  to  the  left  of  zero-lag  on  the  raw  pseudo- 
autocorrelation  function.  This  energy  corresponds  to  the  17  second  peak  found  in  the  original 
spectrum.  The  reason  for  the  misalignment  is  due  to  the  initial  estimate  of  group  velocities 
from  the  MFT.  At  the  18  second  null,  there  is  a  zero  crossing  in  the  complex  spectrum, 
resulting  in  an  instantaneous  phase  change  of  ir  radians.  The  group  delay  at  this  point  should 
be  impulsive,  since  it  is  the  derivative  of  the  phase  with  respect  to  angular  frequency 
(Papoulis,  1962,  p.  134),  and  this  should  be  seen  as  an  impulse  in  group  velocity.  However, 
the  Gaussian  MFT  filter  smooths  the  impulse,  causing  a  distortion  of  the  group  delay  in  the 
adjacent  17  second  peak.  As  a  result,  the  PMF  is  constructed  incorrectly  on  the  first  itera¬ 
tion,  with  the  17  second  energy  to  the  left  of  zero-lag. 

It  is  possible  for  the  PMF  to  correct  the  signal  phase,  but  as  seen  from  Figure  9,  this 
would  require  the  time  window  to  be  at  least  twice  as  wide  as  shown,  reducing  considerably 
the  ability  of  the  filter  to  distinguish  signal  from  noise.  Figure  10  shows  the  effect  of  losing 
the  17  second  peak  on  the  isolated  PMF  and  FVF  modes,  using  the  given  windows. 

Figure  11  shows  the  filtered  spectra  superimposed  on  the  original  spectrum,  with  the 
residuals  between  the  two  plotted  below.  Since  all  three  methods  lost  the  17  second  peak,  the 
residual  plots  start  at  the  18  second  null.  The  MFT  shows  the  effect  of  phase  distortion  on 
the  amplitude  spectrum,  with  the  large  residual  occurring  at  20  seconds.  This  can  be 
confirmed  by  observing  the  near  vertical  slope  to  the  group  velocity  dispersion  in  Figure  8. 
This  again  points  out  the  necessity  for  the  residual  dispersion  method  in  MFT  amplitude 
analysis. 

The  effect  of  curvature  corrections  (19)  on  the  F\T  cosine  windows  can  be  seen  in  Figure 
12.  Maximum  widths  at  18,  20,  and  25  seconds  correspond  to  the  peaks  and  nulls  of  the  spec¬ 
trum  in  Figure  8.  It  would  seem  that  the  largest  width  should  correspond  to  the  point  of 
sharpest  curvature  in  the  spectrum,  at  17  seconds.  However,  a  review  of  equation  (17)  shows 
that  the  bias  estimates  are  constructed  from  the  Parzen  windowed  pseudo-autocorrelation 
function.  Examination  of  Figure  11  again  shows  that  the  window  would  have  to  be  widened 
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considerably  to  extract  this  portion  of  the  spectrum. 

Conclusion 

The  method  of  frequency  variable  filters  (FVF)  is  a  viable  alternative  to  both  the  MFT 
and  PMF  techniques  of  extracting  normal  mode  amplitudes  from  propagating  multi-mode 
surface-waves.  It  has  the  advantage  of  explicitly  addressing  the  problem  of  spectral  bias, 
which  is  an  unavoidable  side  effect  of  any  type  of  convolutional  smoothing  in  the  frequency 
domain.  FVF  will  not  be  successful  on  all  types  of  surface-wave  spectra,  as  indicated  by 
event  III.  This  particular  case  illustrates  the  fundamental  tradeoff  between  frequency  and  time 
domain  resolution.  To  successfully  extract  the  17  second  peak,  the  time  windows  would  have 
to  be  so  wide  that  the  filters  would  be  essentially  allpass.  Both  MFT  and  PMF  had  the  same 
problem  with  this  type  of  event. 

FVF  and  PMF  have  an  advantage  over  MFT,  in  that  the  results  of  processing  can  be 
seen  both  in  the  time  and  frequency  domain.  The  pseudo-autocorrelation  functions  can  be 
used  as  diagnostic  tools,  and  the  final  isolated  mode  can  be  compared  to  the  original  seismo¬ 
gram.  In  addition,  phase  distortion  can  have  a  significant  effect  on  the  MFT,  requiring  appli¬ 
cation  of  a  phase-matched  filter  in  the  form  of  the  residual  dispersion  method  (case  2).  Both 
PMF  and  FVF  remove  the  phase  as  an  automatic  part  of  processing  For  the  above  reasons, 
it  is  recommended  that  the  FVF,  and  in  some  cases  the  PMF,  be  preferred  over  the  MFT  as 
filters  for  isolating  surface-wave  normal  mode  amplitude  spectra. 


Acknowledgements 

™,e  aUth0rS  afe  grateful  10  Hafidh  A  A  Ghalib  and  Hanan  Al-Khatib  for  bringing  the 
WWSSN  seismograms  used  in  this  paper  to  our  attention.  This  work  was  supported  bv  the 
AFSC  under  Contract  F19628-85-K-0029. 

References 

Der,  Z.  A.  (1986).  Comments  on  the  paper  "Estimation  of  scalar  moments  from  explosion- 
generated  surface  waves"  by  Jeffry  L.  Stevens,  Bull.  Seism.  Soe  Am  76  180f>- 
1824. 


Dziewonski,  A.,  S.  Bloch,  and  M.  Landisman  (1969).  A  technique  for  the  analysis  of  transient 
seismic  signals,  Bull.  Seism.  Soc.  Am.  59,  427-444. 

Dziewonski,  A.,  and  A.  Hales  (1972).  Numerical  analysis  of  dispersed  seismic  waves,  in  B.  A. 

Bolt  (editor),  Methods  of  computational  physics,  vol.  11,  Academic  Press  New 
York,  39-85. 

Dziewonski,  A.,  J.  Mills,  and  S.  Bloch  (1972).  Residual  dispersion  measurement  -  a  new 
method  of  surface-wave  analysis,  Bull.  Seism.  Soc.  Am.  62,  129-139. 

Harris,  F.  J.  (1978).  On  the  use  of  windows  for  harmonic  analysis  with  the  discrete  Fourier 
transform,  Proceedings  of  the  IEEE  66,  51-83. 

Herrin,  E.  and  T.  Goforth  (1977).  Phase-matched  filters:  application  to  the  study  of  Rayleigh 
waves,  Bull.  Seis.  Soc.  Am.  67,  1259-1275. 

Herrmann,  R.  B.  (1973).  Some  aspects  of  band-pass  filtering  of  surface  waves,  Bull  Seism 
Soc.  Am.  63,  663-671. 

Hwang,  H.  J.  and  B.  J.  Mitchell  (1986).  Interstation  surface  wave  analysis  by  frequency- 


■■nHBH  IVIA'  UBMHMHI  \  m  Hi  \  aua  iwracnaMt  warn 


domain  Wiener  deconvolution  and  modal  isolation,  Bull.  Seism.  Soc.  ,4m.  76,  847- 
864. 


Jenkins,  Gwilym  M.,  and  D.  G.  Watts  (1068).  Spectral  Analysis  and  Its  Applications,  Holden 
Day,  San  Francisco. 

Landisman,  M.,  A  Dziewonski,  and  Y.  Sato  (1969).  Recent  improvements  in  the  analysis  of 
surface  wave  observations,  Geophys.  J.  17,  369-403. 

Papoulis,  A.  (1962).  The  Fourier  Integral  and  its  Applications,  McGraw-Hill,  New  York. 

Russell,  D  R.,  H.  J.  Hwang,  and  R.  B.  Herrmann  (1986).  Matched  filters  vs.  time  variable 
filters,  Abstract  for  the  AFGL/DARPA  Review  of  Nuclear  Test  Monitoring 
Research,  U.S.  Air  Force  Academy,  Colorado  Springs,  Colorado,  May  6-8  1986. 

Stevens,  J.  L.  (1986a).  Estimation  of  scalar  moments  from  explosion-generated  surface  waves. 
Bull.  Seism.  Soc.  Am.  76,  123-151. 

Stevens,  J.  L.  (1986b).  Reply  to  Z.  Der’s  "Comments  on  the  Paper  'Estimation  of  Scalar 
Moments  from  Explosion-Generated  Surface  Waves’,"  Bull.  Seisin.  Soc.  ,4m.  76. 
1825-1829. 

Yacoub,  N.  K.  (1983).  "Instantaneous  Amplitudes":  a  new  method  to  measure  seismic  magni¬ 
tude,  Bull.  Seism.  Soc.  Am.  7 S,  1345-1355. 


Department  of  Earth  and  Atmospheric  Sciences 

Saint  Louis  University 

P.O.  Box  8099,  Laclede  Station 

St.  Louis,  Missouri  63156 


^  S-  «*•-  M-JV.  A  AV ,V.  /■  /-  -C- 


V*  ■  w  ^  Tl IL  k 


TIME  (*«c) 


FIG.  3.  Comparison  of  the  amplitude  spectra  for  the  isolated  fundamental 
modes.  HEAVY  LINE:  filtered  fundamental  modes  LIGHT  LINE:  theoretical 
fundamental  mode. 
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FIG.  5.  Event  II  pseudo-autocorrelation  functions,  a):  raw  pseudo 
autocorrelation  function,  b):  final  pseudoautocorrelation  function  for  PMF.  c): 
final  pseudoautocorrelation  function  for  FVF.  The  time  marks  indicate  the  two 
sided  window  width. 
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FIG  6.  Event  II  FVF  seismograms,  a):  raw  seismogram,  b):  fundamental 
mode  isolated  by  FVF.  c):  residual  obtained  by  subtracting  the  fundamental 
mode  from  the  original  seismogram. 
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FIG.  8.  Event  III  spectrum  (top)  and  MFT  contour  plot  (bottom).  Time 
domain  seismograms  are  plotted  to  the  right.  Notice  the  non-linear  time  scale  for 
the  MFT  plot.  The  arrow  in  the  upper  left  corner  points  out  the  location  of  the 
17  second  spectra!  peak  referred  to  in  the  text. 
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FIG.  9.  Event  111  pseudo-autocorrelation  functions,  a):  raw  pseudo- 
autocorrelation  function,  b):  final  pseudo-autocorrelation  function  for  PMF.  c). 
final  pseudo-autocorrelation  function  for  FVF.  The  time  marks  indicate  the  two- 
sided  window  width. 
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Fie  11.  Fvcnt  III  amplitude  spectra.  TOP:  comparison  of  the  original 
amplitude  spectrum  with  the  filtered  fundamental  modes.  LIGHT  lines 
correspond  to  the  original  amplitude  spectrum,  and  IILAVA  lines  arc  the  filtered 
modes.  BOTTOM:  residuals  between  the  original  spectrum  and  the  filtered 
modes.  Nolire  that  the  vertical  amplitude  scale  on  these  plots  is  linenr. 
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FIG  1'2.  FVF  and  PMF  one-side  cosine  window  widths  calculated  for  events 
I,  H,  and  III,  respectively.  HEAVY  line.  FVF  wiidow  width.  LIGHT  line:  PMF 
window  width. 
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